text:
!********************************************************************************
!   FPML_QUAD: Fourth order Parallelizable Modification of Laguerre's method in Quadruple Precision
!   Author: Thomas R. Cameron, Davidson College
!   Last Modified: 1 November 2018
!********************************************************************************
module fpml_quad
    use eft, only : dp
    implicit none
    integer, parameter                  :: qp = selected_real_kind(2*precision(1.0_dp))
    real(kind=qp), parameter            :: eps = epsilon(1.0_qp), big = huge(1.0_qp), small = tiny(1.0_qp)

contains
    !************************************************
    !                       main quad               *
    !************************************************
    ! Computes the roots of a polynomial of degree
    ! deg whose coefficients are stored in p. The
    ! root approximations are stored in roots, the
    ! backward error in each approximation in berr,
    ! and the condition number of each root
    ! approximation is stored in cond.
    !************************************************
    subroutine main_quad(poly, deg, roots, berr, cond, conv, itmax)
        implicit none
        ! argument variables
        integer, intent(in)             :: deg, itmax
        integer, intent(out)            :: conv(:)
        real(kind=qp), intent(out)      :: berr(:), cond(:)
        complex(kind=qp), intent(in)    :: poly(:)
        complex(kind=qp), intent(out)   :: roots(:)
        ! local variables
        integer                         :: i, j, nz
        real(kind=qp)                   :: r
        real(kind=qp), dimension(deg+1) :: alpha
        complex(kind=qp)                :: b, c, z
        ! intrinsic functions
        intrinsic                       :: abs, sqrt

        ! precheck
        alpha = abs(poly)
        if(alpha(deg+1)<small) then
            write(*,'(A)') 'Warning: leading coefficient too small.'
            return
        elseif(deg==1) then
            roots(1) = -poly(1)/poly(2)
            conv = 1
            berr = 0.0_qp
            cond(1) = (alpha(1) + alpha(2)*abs(roots(1)))/(abs(roots(1))*alpha(2))
            return
        elseif(deg==2) then
            b = -poly(2)/(2*poly(3))
            c = sqrt(poly(2)**2-4*poly(3)*poly(1))/(2*poly(3))
            roots(1) = b - c
            roots(2) = b + c
            conv = 1
            berr = 0.0_qp
            cond(1) = (alpha(1)+alpha(2)*abs(roots(1))+alpha(3)*abs(roots(1))**2)/(abs(roots(1))* &
                        abs(poly(2)+2.0_qp*poly(3)*roots(1)))
            cond(2) = (alpha(1)+alpha(2)*abs(roots(2))+alpha(3)*abs(roots(2))**2)/(abs(roots(2))* &
                        abs(poly(2)+2.0_qp*poly(3)*roots(2)))
            return
        end if
        ! initial estimates
        conv = (/ (0, i=1,deg)/)
        nz = 0
        call estimates(alpha, deg, roots, conv, nz)
        ! main loop
        alpha = (/ (alpha(i)*(3.8_qp*(i-1)+1),i=1,deg+1)/)
        do i=1,itmax
            do j=1,deg
                if(conv(j)==0) then
                    z = roots(j)
                    r = abs(z)
                    if(r > 1.0_qp) then
                        call rcheck_lag(poly, alpha, deg, b, c, z, r, conv(j), berr(j), cond(j))
                    else
                        call check_lag(poly, alpha, deg, b, c, z, r, conv(j), berr(j), cond(j))
                    end if
                    if(conv(j)==0) then
                        call modify_lag(deg, b, c, z, j, roots)
                        roots(j) = roots(j) - c
                    else
                        nz = nz + 1
                        if(nz==deg) go to 10
                    end if
                end if
            end do
        end do
        ! final check
        10 continue
        if(minval(conv)==1) then
            return
        else
            ! display warrning
            write(*,'(A)') 'Some root approximations did not converge or experienced overflow/underflow.'
            ! compute backward error and condition number for roots that did not converge;
            ! note that this may produce overflow/underflow.
            do j=1,deg
                if(conv(j) .ne. 1) then
                    z = roots(j)
                    r = abs(z)
                    if(r>1.0_qp) then
                        z = 1/z
                        r = 1/r
                        c = 0
                        b = poly(1)
                        berr(j) = alpha(1)
                        do i=2,deg+1
                            c = z*c + b
                            b = z*b + poly(i)
                            berr(j) = r*berr(j) + alpha(i)
                        end do
                        cond(j) = berr(j)/abs(deg*b-z*c)
                        berr(j) = abs(b)/berr(j)
                    else
                        c = 0
                        b = poly(deg+1)
                        berr(j) = alpha(deg+1)
                        do i=deg,1,-1
                            c = z*c + b
                            b = z*b + poly(i)
                            berr(j) = r*berr(j) + alpha(i)
                        end do
                        cond(j) = berr(j)/(r*abs(c))
                        berr(j) = abs(b)/berr(j)
                    end if
                end if
            end do
        end if
    end subroutine main_quad
    !************************************************
    !                       rcheck_lag              *
    !************************************************
    ! Computes backward error of root approximation
    ! with moduli greater than 1.
    ! If the backward error is less than eps, then
    ! both backward error and condition number are
    ! computed. Otherwise, the Laguerre correction terms
    ! are computed and stored in variables b and c.
    !************************************************
    subroutine rcheck_lag(p, alpha, deg, b, c, z, r, conv, berr, cond)
        implicit none
        ! argument variables
        integer, intent(in)             :: deg
        integer, intent(out)            :: conv
        real(kind=qp), intent(in)       :: alpha(:), r
        real(kind=qp), intent(out)      :: berr, cond
        complex(kind=qp), intent(in)    :: p(:), z
        complex(kind=qp), intent(out)   :: b, c
        ! local variables
        integer                         :: k
        real(kind=qp)                   :: rr
        complex(kind=qp)                :: a, zz
        ! intrinsic functions
        intrinsic                       :: abs

        ! evaluate polynomial and derivatives
        zz = 1/z
        rr = 1/r
        a = p(1)
        b = 0
        c = 0
        berr = alpha(1)
        do k=2,deg+1
            c = zz*c + b
            b = zz*b + a
            a = zz*a + p(k)
            berr = rr*berr + alpha(k)
        end do
        ! laguerre correction/ backward error and condition
        if(abs(a)>eps*berr) then
            b = b/a
            c = 2*(c/a)
            c = zz**2*(deg-2*zz*b+zz**2*(b**2-c))
            b = zz*(deg-zz*b)
            if(check_nan_inf(b) .or. check_nan_inf(c)) conv = -1
        else
            cond = berr/abs(deg*a-zz*b)
            berr = abs(a)/berr
            conv = 1
        end if
    end subroutine rcheck_lag
    !************************************************
    !                       check_lag               *
    !************************************************
    ! Computes backward error of root approximation
    ! with moduli less than or equal to 1.
    ! If the backward error is less than eps, then
    ! both backward error and condition number are
    ! computed. Otherwise, the Laguerre correction terms
    ! Gj and Hj are computed and stored in variables
    ! b and c, respectively.
    !************************************************
    subroutine check_lag(p, alpha, deg, b, c, z, r, conv, berr, cond)
        implicit none
        ! argument variables
        integer, intent(in)             :: deg
        integer, intent(out)            :: conv
        real(kind=qp), intent(in)       :: alpha(:), r
        real(kind=qp), intent(out)      :: berr, cond
        complex(kind=qp), intent(in)    :: p(:), z
        complex(kind=qp), intent(out)   :: b, c
        ! local variables
        integer                         :: k
        complex(kind=qp)                :: a
        ! intrinsic functions
        intrinsic                       :: abs

        ! evaluate polynomial and derivatives
        a = p(deg+1)
        b = 0
        c = 0
        berr = alpha(deg+1)
        do k=deg,1,-1
            c = z*c + b
            b = z*b + a
            a = z*a + p(k)
            berr = r*berr + alpha(k)
        end do
        ! laguerre correction/ backward error and condition
        if(abs(a)>eps*berr) then
            b = b/a
            c = b**2 - 2*(c/a)
            if(check_nan_inf(b) .or. check_nan_inf(c)) conv = -1
        else
            cond = berr/(r*abs(b))
            berr = abs(a)/berr
            conv = 1
        end if
    end subroutine check_lag
    !************************************************
    !                       modify_lag              *
    !************************************************
    ! Computes modified Laguerre correction term of
    ! the jth rooot approximation.
    ! The coefficients of the polynomial of degree
    ! deg are stored in p, all root approximations
    ! are stored in roots. The values b, and c come
    ! from rcheck_lag or check_lag, c will be used
    ! to return the correction term.
    !************************************************
    subroutine modify_lag(deg, b, c, z, j, roots)
        implicit none
        ! argument variables
        integer, intent(in)             :: deg, j
        complex(kind=qp), intent(in)    :: roots(:), z
        complex(kind=qp), intent(inout) :: b, c
        ! local variables
        integer                         :: k
        complex(kind=qp)                :: t
        ! intrinsic functions
        intrinsic                       :: abs, sqrt

        ! Aberth correction terms
        do k=1,j-1
            t = 1/(z - roots(k))
            b = b - t
            c = c - t**2
        end do
        do k=j+1,deg
            t = 1/(z - roots(k))
            b = b - t
            c = c - t**2
        end do
        ! Laguerre correction
        t = sqrt((deg-1)*(deg*c-b**2))
        c = b + t
        b = b - t
        if(abs(b)>abs(c)) then
            c = deg/b
        else
            c = deg/c
        end if
    end subroutine modify_lag
    !************************************************
    !                       estimates               *
    !************************************************
    ! Computes initial estimates for the roots of an
    ! univariate polynomial of degree deg, whose
    ! coefficients moduli are stored in alpha. The
    ! estimates are returned in the array roots.
    ! The computation is performed as follows: First
    ! the set (i,log(alpha(i))) is formed and the
    ! upper envelope of the convex hull of this set
    ! is computed, its indices are returned in the
    ! array h (in descending order). For i=c-1,1,-1
    ! there are h(i) - h(i+1) zeros placed on a
    ! circle of radius alpha(h(i+1))/alpha(h(i))
    ! raised to the 1/(h(i)-h(i+1)) power.
    !************************************************
    subroutine estimates(alpha, deg, roots, conv, nz)
        implicit none
        ! argument variables
        integer, intent(in)             :: deg
        integer, intent(inout)          :: conv(:), nz
        real(kind=qp), intent(in)       :: alpha(:)
        complex(kind=qp), intent(inout) :: roots(:)
        ! local variables
        integer                         :: c, i, j, k, nzeros
        real(kind=qp)                   :: a1, a2, ang, r, th
        integer, dimension(deg+1)       :: h
        real(kind=qp), dimension(deg+1) :: a
        real(kind=qp), parameter        :: pi2 = 6.2831853071795865_qp, sigma = 0.7_qp
        ! intrinsic functions
        intrinsic                       :: log, cos, sin, cmplx

        ! Log of absolute value of coefficients
        do i=1,deg+1
            if(alpha(i)>0) then
                a(i) = log(alpha(i))
            else
                a(i) = -1E+30_qp
            end if
        end do
        call conv_hull(deg+1, a, h, c)
        k=0
        th=pi2/deg
        ! Initial Estiamtes
        do i=c-1,1,-1
            nzeros = h(i)-h(i+1)
            a1 = alpha(h(i+1))**(1.0_qp/nzeros)
            a2 = alpha(h(i))**(1.0_qp/nzeros)
            if(a1 .le. a2*small) then
                ! r is too small
                r = 0.0_qp
                nz = nz + nzeros
                conv(k+1:k+nzeros) = -1
                roots(k+1:k+nzeros) = cmplx(0,0,kind=qp)
            else if(a1 .ge. a2*big) then
                ! r is too big
                r = big
                nz = nz+nzeros
                conv(k+1:k+nzeros) = -1
                ang = pi2/nzeros
                do j=1,nzeros
                    roots(k+j) = r*cmplx(cos(ang*j+th*h(i)+sigma),sin(ang*j+th*h(i)+sigma),kind=qp)
                end do
            else
                ! r is just right
                r = a1/a2
                ang = pi2/nzeros
                do j=1,nzeros
                    roots(k+j) = r*cmplx(cos(ang*j+th*h(i)+sigma),sin(ang*j+th*h(i)+sigma),kind=qp)
                end do
            end if
            k = k+nzeros
        end do
    end subroutine estimates
    !************************************************
    !                       conv_hull               *
    !************************************************
    ! Computex upper envelope of the convex hull of
    ! the points in the array a, which has size n.
    ! The number of vertices in the hull is equal to
    ! c, and they are returned in the first c entries
    ! of the array h.
    ! The computation follows Andrew's monotone chain
    ! algorithm: Each consecutive three pairs are
    ! tested via cross to determine if they form
    ! a clockwise angle, if so that current point
    ! is rejected from the returned set.
    !************************************************
    subroutine conv_hull(n, a, h, c)
        implicit none
        ! argument variables
        integer, intent(in)         :: n
        integer, intent(inout)      :: c
        integer, intent(inout)      :: h(:)
        real(kind=qp), intent(in)   :: a(:)
        ! local variables
        integer                     :: i

        ! covex hull
        c=0
        do i=n,1,-1
            do while(c>=2 .and. cross(h, a, c, i)<eps)
                c = c - 1
            end do
            c = c + 1
            h(c) = i
        end do
    end subroutine conv_hull
    !************************************************
    !                       cross                   *
    !************************************************
    ! Returns 2D cross product of OA and OB vectors,
    ! where
    ! O=(h(c-1),a(h(c-1))),
    ! A=(h(c),a(h(c))),
    ! B=(i,a(i)).
    ! If det>0, then OAB makes counter-clockwise turn.
    !************************************************
    function cross(h, a, c, i) result(det)
        implicit none
        ! argument variables
        integer, intent(in)         :: c, i
        integer, intent(in)         :: h(:)
        real(kind=qp), intent(in)   :: a(:)
        ! local variables
        real(kind=qp)               :: det

        ! determinant
        det = (a(i)-a(h(c-1)))*(h(c)-h(c-1)) - (a(h(c))-a(h(c-1)))*(i-h(c-1))
        return
    end function cross
    !************************************************
    !                       check_nan_inf           *
    !************************************************
    ! Check if real or imaginary part of complex
    ! number a is either NaN or Inf.
    !************************************************
    function check_nan_inf(a) result(res)
        implicit none
        ! argument variables
        complex(kind=qp)            :: a
        ! local variables
        logical                     :: res
        real(kind=qp)               :: re_a, im_a
        ! intrinsic functions
        intrinsic                   :: abs

        ! check for nan and inf
        re_a = real(a,kind=qp)
        im_a = aimag(a)
        res = isnan(re_a) .or. isnan(im_a) .or. (abs(re_a)>big) .or. (abs(im_a)>big)
        return
    end function check_nan_inf
end module fpml_quad

max_stars_repo_name: trcameron/FPML-Comp

max_issues_repo_name: trcameron/FPML-Comp

max_forks_repo_name: trcameron/FPML-Comp


submodule(euler_interface_m) euler_prob_0059_m
    implicit none

contains

    module character(len=20) function euler0059()
        write (euler0059, "(i20)") answer()
    end function euler0059

    function answer() result(ret)
        use euler_data_m, only: get_euler_data_0059
        implicit none

        integer :: ret
        integer, allocatable :: encrypted(:), decrypted(:)
        integer, parameter :: n = 26, k = 3
        integer :: letters(n), idx(k), idx2(2*k, k), key(k), i, j
        logical :: next_permutation_avail
        integer :: large_number

        next_permutation_avail = .true.

        call get_euler_data_0059(encrypted)
        large_number = size(encrypted)

        letters = [(i, i=97, 122)]
        allocate (decrypted(size(encrypted)))
        idx = [1, 2, 3]
        idx2(1, :) = [1, 2, 3]
        idx2(2, :) = [1, 3, 2]
        idx2(3, :) = [2, 1, 3]
        idx2(4, :) = [2, 3, 1]
        idx2(5, :) = [3, 1, 2]
        idx2(6, :) = [3, 2, 1]

        ret = 0
        outer: do while (next_permutation_avail)
            inner: do i = 1, size(idx2(:, 1))
                innermost: do j = 1, 3
                    key(j) = letters(idx(idx2(i, j)))
                end do innermost
                call decrypt(encrypted, key, decrypted)
                if (is_english(decrypted)) then
                    ret = sum(decrypted)
                    return
                end if
            end do inner
            call next_permutation(k, n, idx, next_permutation_avail)
        end do outer
    end function answer

    subroutine decrypt(encrypted, key, decrypted)
        integer, intent(in) :: encrypted(:), key(:)
        integer, intent(out) :: decrypted(:)
        integer :: i, k

        do i = 1, size(encrypted)
            k = mod(i, 3)
            if (k == 0) k = 3
            decrypted(i) = ieor(encrypted(i), key(k))
        end do
    end subroutine decrypt

    function is_english(decrypted) result(ret)
        integer, intent(in) :: decrypted(:)
        logical :: ret
        character(len=1) :: check_(8)
        integer :: i, knt(8)

        check_ = ["e", "t", "a", "o", "i", "n", "s", "h"]
        do i = 1, size(check_)
            knt(i) = count(decrypted == iachar(check_(i)))
        end do

        ret = .false.
        if (sum(knt)/real(size(decrypted))*100. > 45. .and. &
            maxloc(knt, dim=1) == 1) ret = .true.
    end function is_english

end submodule euler_prob_0059_m


text:
z <= '1';
x <= a or b;
t <= a or b or c;
y <= a or b;
r <= c or b or a;

ext: f


text:
subroutine mkpwse(igg0,ikk,mode,lossfn, &
& pwse_mass, &
& vol,pi,nwpt,wmax,nbcore,nbocc,ncband,ngkpt,natom,xred,projwf,nlmn, &
& test_bands_se, &
& ntypepaw, &
& pwmatel,tpwmatel, pwjmatel, tpwjmatel, &
& kg,enrgy,cg,npwt,bantot,ncg, &
& indxkpw,indxkbnd,indxkcg,npwarr, &
& kpt,nkpt,nqpt,nsymk,symk,nsym,symrel,syminv, &
& ihlf,lvtrans,bmet,blat,ipaw,itetrahedron, &
& ipwx,ipwndx,npwndx,ntpwndx,napwndx, &
& npwc,npwx,invpw2ndx,pwsymndx,iqsymndx, &
& igmx,igmn,igndx,ikndx,iqndx,isymndx,isymg,npw, &
& nband,nsppol,shiftk,zz,dcse,cse,xse)
implicit none
integer :: mode
! mode = 0 for loss rate
! mode = 1 for power loss
integer :: ipwv,ikk(3),nwpt,nbcore,nbocc,ncband,ngkpt(3),natom,nlmn
integer :: igmn(3),igmx(3)
integer :: npwt,bantot,ncg,nkpt,nqpt,nsym,nsppol
integer :: npw,npwc,npwx,ipw1,npwndx,ntpwndx,napwndx,ipaw,itetrahedron
double precision :: vol,pi,wmax,xred(3,natom)
double precision :: pwse_mass
double complex :: lossfn(nwpt,nqpt+9,ntpwndx)
integer :: test_bands_se(2)
double complex :: projwf(natom,nlmn,nkpt,ncband)
integer :: kg(3,npwt)
double precision :: enrgy(bantot)
double complex :: cg(ncg)
integer :: indxkpw(nkpt),indxkbnd(nkpt)
integer :: indxkcg(nkpt),npwarr(nkpt)
double precision :: kpt(3,nkpt),shiftk(3)
integer :: nsymk(nkpt),symk(nkpt,nsym*2)
integer :: symrel(3,3,nsym),syminv(3,3,nsym)
integer :: lvtrans(3,ngkpt(1),ngkpt(2),ngkpt(3))
integer :: ihlf(nkpt)
double precision :: bmet(3,3),blat(3,3),bmetinv(3,3),bmet_t(3,3),bmet2(3,3)
double precision :: volelmnt, volred
integer :: ipwx(3,npwx),ipwndx(2,napwndx)
integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
integer :: ikndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: isymndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: isymg(3,nkpt,nsym),igsymk(3),igsymq(3)
integer :: nband(nkpt*nsppol)
double complex :: sei,seipw,seiold,seibc,csepw,cseold,csebc
double complex, allocatable :: ssi(:)
double complex :: ssc(-nwpt:nwpt),cse,dcse,cse0,dcse0,cse1(-nwpt:nwpt),zz,ctest,cse2(-nwpt:nwpt)
double precision :: xse,xse0,ssx,ssx1,ssx2,xse2,ssx3
integer :: ii,jj,kk,ll,ix,iy,iz,iskip,ocsign,iiq
integer :: iqq(3),jka(3),jkb(3),jkk(3),iqv(3),ictr(3),iqr(3),ixr(3)
double precision :: qr(3,4),qqr(3,4)
integer :: ixx(3),ixxp(3),ixp2(3),ixp1(3),ixm1(3),ixm2(3)
integer :: ikpt,ikptq
integer :: igg(3),igh(3),igg0(3),igg1(3),igg2(3),idgg1(3),idgg2(3),isgg1(3),isgg(2),igq
integer :: isym,isymq,iqpt,iqsym,iw,iww,ie,je,ies,iqctr,ibp
double precision :: xck(3),xckq(3),qadj(6)
double complex :: cmatel,cmatel2,amatel,amatel2
double complex :: jmatel(3),jmatel2(3),vmatel(3),vmatel2(3)
double complex :: vqmat2(ngkpt(1),ngkpt(2),ngkpt(3))
double complex :: xmat2(ngkpt(1),ngkpt(2),ngkpt(3))
double complex :: jmat2(3,3), rjmat2(3,3), rjmatdum
double precision :: smat2
double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt),vjfactor(3,3,nwpt),vsfactor(3,3,nwpt)
double complex :: vfv(8,nwpt),vfx(8),vsm1(3,3,nwpt),vsm2(3,3,nwpt)
double precision :: vq2(8),vqvtx0(4),vqvtx(4),qkcvt(3,3)
double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3))
double precision :: whi,wlo,wwhi,wwlo,ww,www,wwwlo,wwwhi,enval(8),dw,eshift,evtx0(4),evtx(4)
double precision :: rrpyr(3,4),kvtx(3,4),xk(3,4),ckvtx(3,4),cxk(3,4),rg(3,3),tvol,xpyr0(4),xpyr(4)
double precision :: xkc(3,4),rgc(3,3),tvolc,kvtxc(3,4)
double complex :: vpyr0(4,nwpt),vpyr(4)
double precision :: de21,de31,de32,de41,de42,de43,thresh
double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3)
double complex :: aa0(nwpt),av(3,nwpt)
double precision :: xme,xv(3)
double precision :: avec(3),bgrad(3),xmult
double precision :: abr,rlr,rlr0
integer :: iwh,iwl,ibmin,ibmax,iwhi,iwlo,iehi,ielo
integer :: indxe(4),iwwlim(4)
double precision :: vq(ngkpt(1),ngkpt(2),ngkpt(3)),kf,qq(3),qq2,qp(3),qp2,kv(3),kv2,kmq(3),kmq2,qs(3),ek,ekmq
double precision :: fratio
double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3)),temparray(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: iipw,jjpw,jjpwt
integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx)
integer :: pwsymndx(npwc,2*nsym)
integer :: ipw2,jpw1,jpw2,jw,iqcentr,jqcentr
double precision :: eps1,eps2,eps
double precision :: eval,brd,esprd
double precision :: gfo,gamma(3),pln(3),dist(3)
double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3),stens(3,3)
double complex :: cint
double complex :: wint,wint0,w2int,w2int0,gterm,vcentr
double precision :: wcentr(-nwpt:nwpt),wgrid(nwpt)
integer :: ntypepaw
double complex :: pwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
&                tpwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
&                 pwjmatel(3,ntypepaw,nlmn,nlmn),               &
&                tpwjmatel(3,ntypepaw,nlmn,nlmn)
logical :: lqsing,lqcentr
logical :: lx,lc
double precision :: linterp
double precision :: rr(3,8)
integer :: ivndx(4,6),itet,iv
integer :: nsing
double precision :: wsing(20)
character*4 :: label
external linterp
data rr(1:3,1) /0.0d0,0.0d0,0.0d0/
data rr(1:3,2) /1.0d0,0.0d0,0.0d0/
data rr(1:3,3) /0.0d0,1.0d0,0.0d0/
data rr(1:3,4) /1.0d0,1.0d0,0.0d0/
data rr(1:3,5) /0.0d0,0.0d0,1.0d0/
data rr(1:3,6) /1.0d0,0.0d0,1.0d0/
data rr(1:3,7) /0.0d0,1.0d0,1.0d0/
data rr(1:3,8) /1.0d0,1.0d0,1.0d0/
data ivndx(1:4,1) /1,2,3,5/
data ivndx(1:4,2) /3,5,6,7/
data ivndx(1:4,3) /2,3,5,6/
data ivndx(1:4,4) /2,3,4,6/
data ivndx(1:4,5) /3,4,6,7/
data ivndx(1:4,6) /4,6,7,8/
double precision :: aa,bb,cc,dd,ee,ff
integer :: iter
data iter /30/
logical centerint
integer vcenterint
integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, &
& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, &
& k_ref_DEBUG, i_ref_DEBUG
common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, &
& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, &
& k_ref_DEBUG, i_ref_DEBUG
integer :: idum
double precision :: xdum,xdum2,xdum3,vdum(3),mtxdum(3,3),mtxdum2(3,3)
double complex :: cdum,cdum2,cdum3,cvecdum(3),cvecdum2(3),ctensdum(3,3),ctensdum2(3,3)

!write(6,*) "AAA1"
!xdum = time()

i_prt_DEBUG = 0
j_prt_DEBUG = 0
k_prt_DEBUG = 0
i_ct_DEBUG = 0
j_ct_DEBUG = 0
k_ct_DEBUG = 0
l_ct_DEBUG = 0
i_ref_DEBUG = -1
k_ref_DEBUG = -1
abr=1.d-100
rlr0=3.d-2
rlr=1.d-3
sei=0.d0
ctest=(0.d0,0.d0)
cse=(0.d0,0.d0)
dcse=(0.d0,0.d0)
zz=(0.d0,0.d0)
do ie=-nwpt,nwpt
  cse1(ie)=(0.d0,0.d0)
enddo
xse=0.d0
dw=wmax/dble(nwpt)
!DEBUG
idum = (9.879021/27.2114)/dw+10
!END DEBUG
do iw=1,nwpt
  wgrid(iw)=iw*dw
enddo
volelmnt=(vol*ngkpt(1)*ngkpt(2)*ngkpt(3))/((2.d0*pi)**3)
volred = 1.d0
ikpt=ikndx(ikk(1),ikk(2),ikk(3))
isym=isymndx(ikk(1),ikk(2),ikk(3))
igsymk=isymg(:,ikpt,isym)
do ii=1,3
do jj=1,3
  qkcvt(ii,jj)=blat(ii,jj)/dble(ngkpt(ii))
enddo
enddo
ictr=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/)
iqctr=iqndx(ictr(1),ictr(2),ictr(3))
nsing=0
do ii=1,3
  kv(ii) = dble(ikk(ii)-ictr(ii))/dble(ngkpt(ii)) + igg0(ii)
enddo
kv2 = 0.d0
do ii=1,3
do jj=1,3
  kv2=kv2+kv(ii)*bmet(ii,jj)*kv(jj)
enddo
enddo
ek=kv2/(2.d0*pwse_mass)
kf = (3*pi*pi*(nbocc-nbcore)/vol)**(1.d0/3.d0)
ocsign=-1

do ii=1,3
do kk=1,3
  bmet2(ii,kk) = 0.d0                   ! bmet2 = bmet*bmet
  do jj=1,3
    bmet2(ii,kk) = bmet2(ii,kk) + bmet(ii,jj)*bmet(jj,kk)
  enddo
enddo
enddo

do ix=1,ngkpt(1)
do iy=1,ngkpt(2)
do iz=1,ngkpt(3)
  temparray(ix,iy,iz)=(0.d0,0.d0)
enddo
enddo
enddo

gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0
do ii=1,3
  jj=mod(ii,3)+1
  kk=mod(ii+1,3)+1
  pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2)
  pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1)
  pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1)
  dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln))
enddo
gfo=minval(dist)

!ipaw=0
!write(6,'(a)') '      finding contibutions from:'
!write(6,'(a)') 'band  plane waves    correlation                exchange'
!write(6,'(2f10.6,3x,3f6.2)') ek,9.879021/27.2114,kv
  seibc=sei
  cse2=cse1
  xse2=xse
  do ipw1=1,npwc
!write(6,*) ipw1,npwc
!write(33,*)
!write(33,*) ipw1,npwc
    seipw=sei
    do ie=-nwpt,nwpt
      ssc(ie)=(0.d0,0.d0)
    enddo
    igg1=ipwx(:,ipw1)
    lx=.false.
    lc=ipw1.le.npwc  ! do correlation
    if ((.not.lx).and.(.not.lc)) cycle

! Step 1: Compute matrix elements
    do ix=1,ngkpt(1)
    do iy=1,ngkpt(2)
    do iz=1,ngkpt(3)
!    do ix=1,1
!    do iy=10,10
!    do iz=10,10
      vqmat2(ix,iy,iz)=(0.d0,0.d0)
      ixx=(/ix,iy,iz/)
      iqq=ixx-ictr
      qq = dble(iqq)/dble(ngkpt)
      jka=ikk-iqq
      jkk=mod(jka,ngkpt)
      do ii=1,3
        if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii)
      enddo
      ikptq=ikndx(jkk(1),jkk(2),jkk(3))
      isymq=isymndx(jkk(1),jkk(2),jkk(3))
      igsymq=isymg(:,ikptq,isymq)
      jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ictr
      igg=nint(dble(jkk-jka)/dble(ngkpt))
!write(6,*) '>>>>> iqq = ',iqq
!write(6,*) '>>>>> ikk = ',ikk
!write(6,*) '>>>>> jkk = ',jkk
!write(6,*) '>>>>> jka = ',jka
!write(6,*) '>>>>> igg = ',igg
!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt)
!write(6,'(a,3f10.5)') 'kptq = ',kpt(:,ikptq)
!read(*,*)
      qq=qq+ipwx(:,ipw1)
      kmq = kv - qq
      qq2=0.d0
      kmq2 = 0.d0
      do ii=1,3
      do jj=1,3
        qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj)
        kmq2=kmq2+kmq(ii)*bmet(ii,jj)*kmq(jj)
      enddo
      enddo
      vq(ix,iy,iz)=4.d0*pi/qq2
      ekmq = kmq2/(2.d0*pwse_mass)
!write(6,'(a,3f10.5)') 'qq = ',qq
!write(6,'(a,3f10.5)') 'qq2 = ',qq2
!write(6,*) 'vq = ',vq(ix,iy,iz)
!write(6,*) 'jka = ',jka
!write(6,*) 'jkk = ',jkk
!write(6,*) 'jkb = ',jkb
!write(6,*) 'igg = ',igg
!write(6,*) 'ikpt = ',ikpt
!write(6,*) 'ikptq = ',ikptq
      omega(ix,iy,iz)=ekmq-ek+dw/2
!write(33,'(3i3,2x,4f10.6)') ixx,ek,ekmq,ek-ekmq,9.879021/27.2114
      lqsing=.false.
      xmat2(ix,iy,iz)=(1.d0,0.d0)
      if (iqq(1).eq.0.and.iqq(2).eq.0.and.iqq(3).eq.0) then
        jmat2 = bmetinv*4.d0*pi
        if (mode.eq.1) jmat2 = jmat2*(-omega(ix,iy,iz))
        smat2=0.d0
        if (ipw1.eq.1) then
          smat2=4.d0*pi
          lqsing=.true.
          if (mode.eq.1) smat2 = smat2*(-omega(ix,iy,iz))
        endif
      else
        vqmat2(ix,iy,iz)=xmat2(ix,iy,iz)*vq(ix,iy,iz)
        if (mode.eq.1) vqmat2(ix,iy,iz) = vqmat2(ix,iy,iz)*(-omega(ix,iy,iz))
      endif
!if (ipw1.eq.59.and.(ix.eq.7.or.ix.eq.8).and.(iy.eq.9.or.iy.eq.10).and.(iz.eq.7.or.iz.eq.8)) then
!endif
    enddo
    enddo
    enddo
!read(*,*)

! DEBUG
!write(6,*) "Matrix elements computed"

    if (itetrahedron.eq.0) then
      write(6,*) "Sum over points not yet implemented in self energy calculation"
      write(6,*) "set itetrahedron=1 and run again"
    else
! tetrahedron integration

! Don't waste time integrating over energies with no contribution
      do ix=1,ngkpt(1)
      do iy=1,ngkpt(2)
      do iz=1,ngkpt(3)
        ixx=(/ix,iy,iz/)
        call fhilo(omega,ixx,ngkpt,whi,wlo)
        if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
          wwhi=whi
          wwlo=wlo
        else
          wwhi=max(wwhi,whi)
          wwlo=min(wwlo,wlo)
        endif
      enddo
      enddo
      enddo

      iwlo=nint(wwlo*dble(nwpt)/wmax)
      iwhi=nint(wwhi*dble(nwpt)/wmax)
      if (ocsign.lt.0.d0) then
        ielo=iwlo+1
        iehi=iwhi+nwpt
      else
        ielo=iwlo-nwpt
        iehi=iwhi-1
      endif
      allocate (ssi(ielo:iehi))
      do ie=ielo,iehi
        ssi(ie)=0.d0
      enddo
      ssx=0.d0

      if (lc) then
        do iw=1,nwpt
          do ix=1,ngkpt(1)
          do iy=1,ngkpt(2)
          do iz=1,ngkpt(3)
            ixx=(/ix,iy,iz/)
            iqpt=iqndx(ixx(1),ixx(2),ixx(3))
            call locateelement(ixx,ipw1,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
            if (iqpt.eq.iqctr) then
              vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0)
            else if (jjpw.ne.0) then
              vfactor(ixx(1),ixx(2),ixx(3),iw)= &
&                 vqmat2(ixx(1),ixx(2),ixx(3))* &
&                 (lossfn(iw,iqpt,jjpw)-conjg(lossfn(iw,iqpt,jjpw)))*0.5d0
            else
              vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0)
            endif
!if (ipw1.eq.59.and.iw.eq.247.and.(ix.eq.7.or.ix.eq.8).and.(iy.eq.9.or.iy.eq.10).and.(iz.eq.7.or.iz.eq.8)) then
!write(34,*) ixx,vqmat2(ix,iy,iz),lossfn(iw,iqpt,jjpw)-conjg(lossfn(iw,iqpt,jjpw)),vfactor(ixx(1),ixx(2),ixx(3),iw)
!endif
!if (iw.eq.idum) write(33,'(3i3,2x,2f10.6,2x,f10.6)') ixx,lossfn(iw,iqpt,jjpw),omega(ix,iy,iz)
          enddo
          enddo
          enddo
!write(6,*) "iw = ", iw
          call locateelement(ictr,ipw1,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
          do ii=1,3
          do jj=1,3
            if (jjpw.ne.0) then
              vjfactor(ii,jj,iw)=(0.d0,0.d0)
              do kk=1,3
              do ll=1,3
                iiq=ii + 3*(kk-1)
                vjfactor(ii,jj,iw)= &
&                     vjfactor(ii,jj,iw) + &
&                     0.5d0*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpw))) &
&                     *bmet(kk,ll)*jmat2(ll,jj)
              enddo
              enddo
              iiq=ii + 3*(jj-1)
              vsfactor(ii,jj,iw)=smat2*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpw)))*0.5d0
            else
              vjfactor(ii,jj,iw)=(0.d0,0.d0)
              vsfactor(ii,jj,iw)=(0.d0,0.d0)
            endif
          enddo
          enddo
!do ii=1,3
!write(6,*) dble(vsfactor(ii,:,iw))
!enddo
!write(6,*)
!do ii=1,3
!write(6,*) dimag(vsfactor(ii,:,iw))
!enddo
!read(*,*)
          do ii=1,3
          do jj=1,3
            vsm1(ii,jj,iw)=(0.d0,0.d0)
            do kk=1,3
              vsm1(ii,jj,iw)=vsm1(ii,jj,iw)+bmet(ii,kk)*vsfactor(kk,jj,iw)
            enddo
          enddo
          enddo
          do ii=1,3
          do jj=1,3
            vsm2(ii,jj,iw)=(0.d0,0.d0)
            do kk=1,3
              vsm2(ii,jj,iw)=vsm2(ii,jj,iw)+vsm1(ii,kk,iw)*bmet(kk,jj)
            enddo
          enddo
          enddo
        enddo
      endif

! DEBUG
!write(6,*) "Beginning integration"

! Do the integration
      do ix=1,ngkpt(1)
      do iy=1,ngkpt(2)
      do iz=1,ngkpt(3)
!xdum = 0.d0;
!xdum2 = 0.d0;
!xdum3 = 0.d0;
        ixx=(/ix,iy,iz/)
        iqq=ixx-ictr
        ssx1=ssx
        call fval(omega,ixx,ixxp,ngkpt,enval)
        if (lqsing) call fval(vq,ixx,ixxp,ngkpt,vq2)
        if (lc) then
          do iw=1,nwpt
            call fpol(vfactor(:,:,:,iw),ngkpt,ixx,ixxp,vfv(:,iw))
          enddo
        endif
        if (lx) call fpol(vqmat2(:,:,:),ngkpt,ixx,ixxp,vfx)
        do itet=1,6
! DEBUG
!cdum = 0.d0
          centerint = .false.
          do iv=1,4
            evtx0(iv)=enval(ivndx(iv,itet))
            rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet))
            if (lc) then
              do iw=1,nwpt
                vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw)
              enddo
            endif
            if (lx) xpyr0(iv)=vfx(ivndx(iv,itet))
            if ((ixx(1)+rrpyr(1,iv).eq.ictr(1)).and.(ixx(2)+rrpyr(2,iv).eq.ictr(2)).and.(ixx(3)+rrpyr(3,iv).eq.ictr(3))) then
              centerint = .true.  ! does tetrahedron contain center vertex?
              vcenterint = iv     ! index of center vertex
            endif
          enddo
!write(6,*) ixx
!write(6,*) itet,centerint
!if (centerint) then
!write(6,*) itet
!do iv=1,4
!write(6,*) ixx+rrpyr(:,iv)
!enddo
!!read(*,*)
!endif
          call indxhpsort(4,4,evtx0,indxe)
          evtx=evtx0(indxe)
          do ii=1,4
            jj=indxe(ii)
            kvtx(:,ii) = (iqq(:) + rrpyr(:,jj))
            do kk=1,3
              kvtxc(kk,ii)=dot_product(iqq+rrpyr(:,jj),qkcvt(:,kk))
            enddo
          enddo
          xk(1:3,1)=(/0.d0,0.d0,0.d0/)
          xkc(1:3,1)=(/0.d0,0.d0,0.d0/)
          do ii=1,3
            xk(:,ii+1)=rrpyr(:,indxe(ii+1))-rrpyr(:,indxe(1))
          enddo
          do ii=2,4
            xkc(1:3,ii)=kvtxc(1:3,ii)-kvtxc(1:3,1)
          enddo
          do jj=1,3  ! make contragredient, reduced coordinates
            call cross(xk(:,mod(jj,3)+2),xk(:,mod(jj+1,3)+2),rg(:,jj))
          enddo
          do jj=1,3  ! make contragredient, Cartesian coordinates
            call cross(xkc(:,mod(jj,3)+2),xkc(:,mod(jj+1,3)+2),rgc(:,jj))
          enddo
          tvol=dot_product(xk(1:3,2),rg(1:3,1))
          rg=rg/tvol
          tvol=abs(tvol)
          tvolc=dot_product(xkc(1:3,2),rgc(1:3,1))
          rgc=rgc/tvolc
          tvolc=abs(tvolc)
          iwwlim=nint(evtx/dw)
!if (centerint) write(6,'(i3,3x,3i3)') itet,ixx
!if (centerint) write(6,*) centerint,lqsing,lx
!if (ek.gt.0.2.and.ipw1.eq.1) write(34,*) "A    ",itet,ixx,ssi(1)
          if (lqsing .and. centerint) then
! singular part
! DEBUG
!xdum = 0.d0
!write(6,*)
!write(6,*) "singular integration"
!write(6,*) 'ixx: ',ixx
!write(6,*) 'itet: ',itet,'  centerint: ',centerint
!write(6,*) 'iwwlim: ',iwwlim
            do iv=1,4
! subtraction of integrated singular function from values at vertices
              ixr=ixx+rrpyr(:,iv)
              iqr = ixr - ictr
              qr(:,iv) = dble(iqr)
              if (ixr(1).eq.ictr(1).and.ixr(2).eq.ictr(2).and.ixr(3).eq.ictr(3)) then
                if (lc) then
                  do iw=1,nwpt
                    vpyr0(iv,iw) = 0.d0
                  enddo
                endif
                if (lx) xpyr0(iv) = 0.d0
              else
                qq = qr(:,iv)/dble(ngkpt)
                qqr(:,iv)=qq
                qq2 = 0.d0
                do ii=1,3
                do jj=1,3
                  qq2 = qq2 + qq(ii)*bmet(ii,jj)*qq(jj)
                enddo
                enddo
                if (lc) then
                  do iw=1,nwpt
                    do ii=1,3
                    do jj=1,3
                      vpyr0(iv,iw)=vpyr0(iv,iw)-qq(ii)*vsm2(ii,jj,iw)*qq(jj)/(qq2**2)
                    enddo
                    enddo
                  enddo
                endif
                if (lx) then
                  xpyr0(iv)=xpyr0(iv)-smat2/qq2
                endif
              endif
            enddo
            bgrad=(/0.d0,0.d0,0.d0/)   ! energy gradient
            do jj=1,3
              bgrad=bgrad+(evtx(jj+1)-evtx(1))*rgc(:,jj)
            enddo
            xmult=1.d0/sqrt(dot_product(bgrad,bgrad))
            if (lc) then
              do iww=iwwlim(1)-1,iwwlim(4)+1
                wwwlo=dw*(iww-0.5d0)
                call singular_adaptive_tetint(qqr,evtx,wwwlo,dw,vcenterint,bmet,abr,rlr,rlr0,stens)
                do iw=1,nwpt
                  ie=iww-ocsign*iw
                  if (ie.ge.ielo.and.ie.le.iehi) then
                    do ii=1,3
                    do jj=1,3
                      ssi(ie)=ssi(ie)+vsm2(ii,jj,iw)*stens(jj,ii)
! DEBUG
!if (ie.eq.0.or.ie.eq.1) xdum3=xdum3+dimag(vsm2(ii,jj,iw)*stens(jj,ii))
                    enddo
                    enddo
                  endif
                enddo
              enddo
            endif
            if (lx) then
! DEBUG
!cdum = 0.d0
              wwwlo = min(evtx(1),evtx(2),evtx(3),evtx(4))
              dw = max(evtx(1),evtx(2),evtx(3),evtx(4))-wwwlo
              call singular_adaptive_tetint(qqr,evtx,wwwlo,dw,vcenterint,bmet,abr,rlr,rlr0,stens)
              do ii=1,3
              do jj=1,3
                ssx=ssx+smat2*bmet(ii,jj)*stens(jj,ii)
!! DEBUG
!cdum = cdum+smat2*bmet(ii,jj)*stens(jj,ii)
              enddo
              enddo
!write(6,*) ssx
!write(6,'(i3,2x,3i3,2x,2f12.8)') itet,ixx,cdum
            endif
          elseif (centerint) then
            if (lc) then
              do iww=iwwlim(1)-1,iwwlim(4)+1
                wwwlo=dw*(iww-0.5d0)
                do iw=1,nwpt
                  call vnitetint(rrpyr,evtx0,wwwlo,dw,vjfactor(:,:,iw),bmet,vcenterint,cint)
                  ie=iww-ocsign*iw
                  if (ie.ge.ielo.and.ie.le.iehi) ssi(ie)=ssi(ie)+cint/volelmnt
! DEBUG
!if (ie.eq.0.or.ie.eq.1) xdum2=xdum2+dimag(cint/volelmnt)
                enddo
              enddo
            endif
            if (lx) then
              call vnitetint(rrpyr,evtx0,evtx0(indxe(1)),evtx0(indxe(4))-evtx0(indxe(1)),jmat2,bmet,vcenterint,cint)
              ssx=ssx+cint/volelmnt
            endif
          endif
!if (ek.gt.0.2.and.ipw1.eq.1) write(34,*) "     ",itet,ixx,ssi(1)
! can use basic tetrahedron integration for everything else
! DEBUG
!cdum = 0.d0
!if ((ix.eq.1.or.ix.eq.2).and.(iy.eq.9.or.iy.eq.10).and.(iz.eq.8.or.iz.eq.9)) then
!write(6,*) "vpyr0", vpyr0(:,idum)
!write(6,*) "evtx",evtx
!write(6,*) "www",idum*dw
!write(6,*) "iw",idum
!write(6,*) "iwwlim",iwwlim
!read(*,*)
!endif
          if (lc) then
            do iw=1,nwpt
              aa0(iw)=vpyr0(indxe(1),iw)   ! base value of function
              av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function
              do jj=1,3
                av(1:3,iw)=av(1:3,iw) &
&                      +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj)
              enddo
            enddo
            do iww=iwwlim(1)-1,iwwlim(4)+1
              wwwlo=dw*(iww-0.5d0)
!              wwwhi=dw*(iww+0.5d0)
              call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,wwwlo,dw,sint1,svec)
              do iw=1,nwpt
                ie=iww-ocsign*iw
                if (ie.ge.ielo.and.ie.le.iehi) ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt
! DEBUG
!if (ipw1.eq.59.and.(ix.eq.7).and.(iy.eq.9).and.(iz.eq.7)) then
!if ((ie.eq.0.or.ie.eq.1).and.iw.eq.247) then
!write(34,*) ixx,iw,wwwlo+dw/2.d0,(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt
!write(34,*) "                  ",vpyr0(indxe(1),iw)
!write(34,*) "                  ",vpyr0(indxe(2),iw)
!write(34,*) "                  ",vpyr0(indxe(3),iw)
!write(34,*) "                  ",vpyr0(indxe(4),iw)
!endif
!endif
!if (ie.eq.0.or.ie.eq.1) xdum = xdum+dimag((aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt)
!if (ie.eq.0.or.ie.eq.1) then
!if (ix.eq.1.and.iy.eq.1.and.iz.eq.4.and.itet.eq.1) then
!if ((ix.eq.1.or.ix.eq.2).and.(iy.eq.9.or.iy.eq.10).and.(iz.eq.8.or.iz.eq.9)) then
!if (iw.eq.idum) then
!write(6,*) "aa0",aa0(iw)
!write(6,*) "av",av(:,iw)
!write(6,*) "sint1",sint1
!write(6,*) "svec",svec
!write(6,*) "evtx",evtx
!write(6,*) "wwwlo",wwwlo
!write(6,*) "wwwhi",wwwlo+dw
!write(6,*) "vpyr0", vpyr0(:,iw)
!write(6,*) "vfv", vfv(:,iw)
!read(*,*)
!endif
!endif
!endif
              enddo
            enddo
          endif
! DEBUG
!write(32,'(4i3,2x,f14.10,2x,f14.10)') ixx,itet,(cdum/2.d0)
          if (lx) then
            xme=xpyr0(indxe(1))
            xv=(/0.d0,0.d0,0.d0/)
            do jj=1,3
              xv=xv+(xpyr0(indxe(jj+1))-xpyr0(indxe(1)))*rg(:,jj)
            enddo
            call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,evtx(1),evtx(4)-evtx(1),sint1,svec)
            ssx=ssx+(xme*sint1+dot_product(svec,xv))/volelmnt
! DEBUG
!write(32,'(4i3,2x,f16.10)') ixx,itet,(xme*sint1+dot_product(svec,xv))/volelmnt
!write(32,'(4i3,2x,f14.10,2x,3f14.10)') ixx,itet, &
!& (xme*sint1+dot_product(svec,xv))/volelmnt, &
!& (xpyr0(indxe(1))+xpyr0(indxe(2))+xpyr0(indxe(3))+xpyr0(indxe(4))) &
!& /(volelmnt*6.d0*4.d0)
          endif
!write(33,'(3i3,2x,i3,2x,4f10.6,2x,2f10.6)') ixx,itet,evtx,cdum
!if (ek.gt.0.2.and.ipw1.eq.1) write(34,*) "     ",itet,ixx,ssi(1)
        enddo
! DEBUG
!write(32,'(3i3,2x,4f16.10)') ixx,ssi(0),ssi(1)
!write(33,*) ixx,xdum,xdum2,xdum3
      enddo
      enddo
      enddo
!    if (lc) ssi=-ssi/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
!    if (lx) ssx=-ssx/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol)
      if (lc) ssi=-ocsign*ssi/((2*pi)**3*pi)
      if (lx) ssx=-ssx/((2*pi)**3)

! DEBUG
!write(6,*) "Integration complete"

      if (lc) then
        do ie=-nwpt,nwpt
          eps1=dble(ie)*dw-dw/2
!        ssc(ie)=(0.d0,0.d0)
          do je=ielo,iehi
            if (ie.eq.je) then
              ssc(ie)=ssc(ie)+pi*ssi(je)
            else
              eps2=dble(je)*dw-dw/2
              ssc(ie)=ssc(ie)-(0.d0,1.d0)*ocsign*ssi(je)*dw/(eps2-eps1)
            endif
!if (ie.eq.0) write(16,'(2i4,2f10.3,es12.3,2(2x,2es11.3))') ie,je,eps1*27.2114,eps2*27.2114,1/(eps2-eps1),ssc(ie)*27.2114,dble(ssi(je)*27.2114)
!if (ie.eq.0) write(16,'(2x,i4,2x,f10.3,4x,f10.6,2x,f10.6)') je,eps2*27.2114,ssi(je)*27.2114
          enddo
!if (ipw1.eq.1) write(34,*) ie, ssi(ie)
        enddo
      endif
!do je=ielo,iehi
!write(6,*) je,ssi(je)
!enddo
      deallocate (ssi)
!write(6,'("lstv",2x,f14.10,4x,"(",f14.10,",",f14.10,")")') ssx,(ssc(0)+ssc(1))/2.

      if (lc) cse1=cse1+ssc
      if (lx) xse=xse+ssx
    endif

! DEBUG
  enddo

if (itetrahedron.eq.1) then
  cse=(cse1(0)+cse1(1))/2
  dcse=(cse1(1)-cse1(0))/dw
endif
zz=1.d0/(1.d0-dcse)
!write(6,*) cse*27.2114
!write(6,*) xse*27.2114
!write(6,*) zz
!write(6,*) dble(zz)*dimag(cse)*27.2114

!temparray=log10(temparray)
!temparray=10**(temparray-int(temparray))
!do ix=1,ngkpt(1)
!  write(76,*) ix,'--------------------------------------------------'
!  do iy=1,ngkpt(2)
!!    write(76,'(10es8.1)') dble(temparray(ix,iy,:))
!    write(76,'(10f8.5)') dble(temparray(ix,iy,:))
!  enddo
!enddo

!write(6,*) "AAA2", time()-xdum

return
end subroutine mkpwse

!*************************************************************************

subroutine mkxtpwse(igg0,ikk,mode,xtlossfn, &
& pwse_mass, &
& vol,pi,nwpt,wmax,nbcore,nbocc,ncband,ngkpt,natom,xred,projwf,nlmn, &
& test_bands_se, &
& ntypepaw, &
& pwmatel,tpwmatel, pwjmatel, tpwjmatel, &
& kg,enrgy,cg,npwt,bantot,ncg, &
& indxkpw,indxkbnd,indxkcg,npwarr, &
& kpt,nkpt,nqpt,nxtqpt,nsymk,symk,nsym,symrel,syminv, &
& ihlf,lvtrans,bmet,blat,ipaw,itetrahedron, &
& ipwx,ipwndx,npwndx,ntpwndx,napwndx, &
& pwlo,pwhi,npwc,npwx,invpw2ndx,pwsymndx,iqsymndx, &
& igmx,igmn,igndx,ikndx,iqndx,isymndx,isymg,npw, &
& nxtgkpt,xtqpt,xtqpt2qpt,qpt2xtqpt,qpt2xtg,ixtqndx, &
& nband,nsppol,shiftk,zz,dcse,cse,xse)
implicit none
integer :: mode
! mode = 0 for loss rate
! mode = 1 for power loss
integer :: ipwv,ikk(3),nwpt,nbcore,nbocc,ncband,ngkpt(3),natom,nlmn
integer :: nxtgkpt(3)
integer :: igmn(3),igmx(3)
integer :: npwt,bantot,ncg,nkpt,nqpt,nxtqpt,nsym,nsppol
integer :: npw,npwc,npwx,pwlo,pwhi,ipw1,npwndx,ntpwndx,napwndx,ipaw,itetrahedron
double precision :: vol,pi,wmax,xred(3,natom)
double precision :: pwse_mass
double complex :: xtlossfn(nwpt,nxtqpt,npwx)
integer :: test_bands_se(2)
double complex :: projwf(natom,nlmn,nkpt,ncband)
integer :: kg(3,npwt)
double precision :: enrgy(bantot)
double complex :: cg(ncg)
integer :: indxkpw(nkpt),indxkbnd(nkpt)
integer :: indxkcg(nkpt),npwarr(nkpt)
double precision :: kpt(3,nkpt),shiftk(3)
integer :: nsymk(nkpt),symk(nkpt,nsym*2)
integer :: symrel(3,3,nsym),syminv(3,3,nsym)
integer :: lvtrans(3,ngkpt(1),ngkpt(2),ngkpt(3))
integer :: ihlf(nkpt)
double precision :: bmet(3,3),blat(3,3),bmetinv(3,3),bmet_t(3,3),bmet2(3,3)
double precision :: volelmnt, volred
integer :: ipwx(3,npwx),ipwndx(2,napwndx)
integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
integer :: ikndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: ixtqndx(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3))
integer :: isymndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: isymg(3,nkpt,nsym),igsymk(3),igsymq(3)
double precision :: xtqpt(3,nxtqpt)
integer :: xtqpt2qpt(nxtqpt),qpt2xtqpt(nqpt),qpt2xtg(nqpt,3)
integer :: nband(nkpt*nsppol)
double complex :: sei,seipw,seiold,seibc,csepw,cseold,csebc
double complex, allocatable :: ssi(:)
double complex :: ssc(-nwpt:nwpt),cse,dcse,cse0,dcse0,cse1(-nwpt:nwpt),zz,ctest,cse2(-nwpt:nwpt)
double precision :: xse,xse0,ssx,ssx1,ssx2,xse2,ssx3
integer :: ii,jj,kk,ll,ix,iy,iz,iskip,ocsign,iiq
integer :: iqq(3),jka(3),jkb(3),jkk(3),iqv(3),ictr(3),ixtctr(3),iqr(3),ixr(3)
double precision :: qr(3,4),qqr(3,4)
integer :: ixx(3),ixt(3),ixxp(3),ixtp(3),ixp2(3),ixp1(3),ixm1(3),ixm2(3)
integer :: ikpt,ikptq
integer :: igg(3),igh(3),igg0(3),igg1(3),igg2(3),idgg1(3),idgg2(3),isgg1(3),isgg(2),igq
integer :: isym,isymq,iqpt,ixtqpt,iqsym,iw,iww,ie,je,ies,iqctr,ibp
double precision :: xck(3),xckq(3),qadj(6)
double complex :: cmatel,cmatel2,amatel,amatel2
double complex :: jmatel(3),jmatel2(3),vmatel(3),vmatel2(3)
double complex :: vqmat2(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3))
double complex :: xmat2(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3))
double complex :: jmat2(3,3), rjmat2(3,3), rjmatdum
double precision :: smat2
double complex :: vfactor(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3),nwpt)
double complex :: vfv(8,nwpt),vfx(8)
double precision :: vq2(8),vqvtx0(4),vqvtx(4),qkcvt(3,3)
double precision :: omega(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3))
double precision :: whi,wlo,wwhi,wwlo,ww,www,wwwlo,wwwhi,enval(8),dw,eshift,evtx0(4),evtx(4)
double precision :: rrpyr(3,4),kvtx(3,4),xk(3,4),ckvtx(3,4),cxk(3,4),rg(3,3),tvol,xpyr0(4),xpyr(4)
double precision :: xkc(3,4),rgc(3,3),tvolc,kvtxc(3,4)
double complex :: vpyr0(4,nwpt),vpyr(4)
double precision :: de21,de31,de32,de41,de42,de43,thresh
double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3)
double complex :: aa0(nwpt),av(3,nwpt)
double precision :: xme,xv(3)
double precision :: avec(3),bgrad(3),xmult
double precision :: abr,rlr,rlr0
integer :: iwh,iwl,ibmin,ibmax,iwhi,iwlo,iehi,ielo
integer :: indxe(4),iwwlim(4)
double precision :: vq(nxtgkpt(1),nxtgkpt(2),nxtgkpt(3))
double precision :: kf,qq(3),qq2,qp(3),qp2,kv(3),kv2,kmq(3),kmq2,qs(3),ek,ekmq
double precision :: fratio
double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3)),temparray(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: iipw,jjpw,jjpwt
integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx)
integer :: pwsymndx(npwc,2*nsym)
integer :: ipw2,jpw1,jpw2,jw,iqcentr,jqcentr
double precision :: eps1,eps2,eps
double precision :: eval,brd,esprd
double precision :: gfo,gamma(3),pln(3),dist(3)
double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3),stens(3,3)
double complex :: cint
double complex :: wint,wint0,w2int,w2int0,gterm,vcentr
double precision :: wcentr(-nwpt:nwpt),wgrid(nwpt)
integer :: ntypepaw
double complex :: pwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
&                tpwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
&                 pwjmatel(3,ntypepaw,nlmn,nlmn),               &
&                tpwjmatel(3,ntypepaw,nlmn,nlmn)
logical :: lqsing,lqcentr
logical :: lx,lc
double precision :: linterp
double precision :: rr(3,8)
integer :: ivndx(4,6),itet,iv
integer :: nsing
double precision :: wsing(20)
character*4 :: label
external linterp
data rr(1:3,1) /0.0d0,0.0d0,0.0d0/
data rr(1:3,2) /1.0d0,0.0d0,0.0d0/
data rr(1:3,3) /0.0d0,1.0d0,0.0d0/
data rr(1:3,4) /1.0d0,1.0d0,0.0d0/
data rr(1:3,5) /0.0d0,0.0d0,1.0d0/
data rr(1:3,6) /1.0d0,0.0d0,1.0d0/
data rr(1:3,7) /0.0d0,1.0d0,1.0d0/
data rr(1:3,8) /1.0d0,1.0d0,1.0d0/
data ivndx(1:4,1) /1,2,3,5/
data ivndx(1:4,2) /3,5,6,7/
data ivndx(1:4,3) /2,3,5,6/
data ivndx(1:4,4) /2,3,4,6/
data ivndx(1:4,5) /3,4,6,7/
data ivndx(1:4,6) /4,6,7,8/
double precision :: aa,bb,cc,dd,ee,ff
integer :: iter
data iter /30/
logical centerint
integer vcenterint
integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, &
& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, &
& k_ref_DEBUG, i_ref_DEBUG
common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, &
& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, &
& k_ref_DEBUG, i_ref_DEBUG
integer :: idum
double precision :: xdum,xdum2,xdum3,vdum(3),mtxdum(3,3),mtxdum2(3,3)
double complex :: cdum,cdum2,cdum3,cvecdum(3),cvecdum2(3),ctensdum(3,3),ctensdum2(3,3)

!write(6,*) "AAA1"
!xdum = time()
!write(6,*) xse

i_prt_DEBUG = 0
j_prt_DEBUG = 0
k_prt_DEBUG = 0
i_ct_DEBUG = 0
j_ct_DEBUG = 0
k_ct_DEBUG = 0
l_ct_DEBUG = 0
i_ref_DEBUG = -1
k_ref_DEBUG = -1
abr=1.d-100
rlr0=3.d-2
rlr=1.d-3
sei=0.d0
ctest=(0.d0,0.d0)
cse=(0.d0,0.d0)
dcse=(0.d0,0.d0)
zz=(0.d0,0.d0)
do ie=-nwpt,nwpt
  cse1(ie)=(0.d0,0.d0)
enddo
xse=0.d0
dw=wmax/dble(nwpt)
!DEBUG
idum = (9.879021/27.2114)/dw+10
!END DEBUG
do iw=1,nwpt
  wgrid(iw)=iw*dw
enddo
volelmnt=(vol*nxtgkpt(1)*nxtgkpt(2)*nxtgkpt(3))/((2.d0*pi)**3)
volred = 1.d0
ikpt=ikndx(ikk(1),ikk(2),ikk(3))
isym=isymndx(ikk(1),ikk(2),ikk(3))
igsymk=isymg(:,ikpt,isym)
do ii=1,3
do jj=1,3
  qkcvt(ii,jj)=blat(ii,jj)/dble(nxtgkpt(ii))
enddo
enddo
ictr=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/)
ixtctr=(/nxtgkpt(1)/2,nxtgkpt(2)/2,nxtgkpt(3)/2/)
iqctr=iqndx(ictr(1),ictr(2),ictr(3))
nsing=0
do ii=1,3
  kv(ii) = dble(ikk(ii)-ictr(ii))/dble(ngkpt(ii)) + igg0(ii)
enddo
kv2 = 0.d0
do ii=1,3
do jj=1,3
  kv2=kv2+kv(ii)*bmet(ii,jj)*kv(jj)
enddo
enddo
ek=kv2/(2.d0*pwse_mass)
kf = (3*pi*pi*(nbocc-nbcore)/vol)**(1.d0/3.d0)
ocsign=-1

do ii=1,3
do kk=1,3
  bmet2(ii,kk) = 0.d0                   ! bmet2 = bmet*bmet
  do jj=1,3
    bmet2(ii,kk) = bmet2(ii,kk) + bmet(ii,jj)*bmet(jj,kk)
  enddo
enddo
enddo

!do ix=1,ngkpt(1)
!do iy=1,ngkpt(2)
!do iz=1,ngkpt(3)
!  temparray(ix,iy,iz)=(0.d0,0.d0)
!enddo
!enddo
!enddo

gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0
do ii=1,3
  jj=mod(ii,3)+1
  kk=mod(ii+1,3)+1
  pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2)
  pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1)
  pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1)
  dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln))
enddo
gfo=minval(dist)

!ipaw=0
!write(6,'(a)') '      finding contibutions from:'
!write(6,'(a)') 'band  plane waves    correlation                exchange'
!write(6,'(2f10.6,3x,3f6.2)') ek,9.879021/27.2114,kv
  seibc=sei
  cse2=cse1
  xse2=xse
  do ipw1=pwlo,pwhi
    seipw=sei
    do ie=-nwpt,nwpt
      ssc(ie)=(0.d0,0.d0)
    enddo
    igg1=ipwx(:,ipw1)
    lx=.false.
    lc=ipw1.le.npwx  ! do correlation
    if ((.not.lx).and.(.not.lc)) cycle

! Step 1: Compute matrix elements
    do ix=1,nxtgkpt(1)
    do iy=1,nxtgkpt(2)
    do iz=1,nxtgkpt(3)
!    do ix=1,1
!    do iy=10,10
!    do iz=10,10
      vqmat2(ix,iy,iz)=(0.d0,0.d0)
      ixt=(/ix,iy,iz/)
      ixx=ixt*ngkpt/nxtgkpt
      iqq=ixx-ictr
      qq = dble(iqq)/dble(ngkpt)
      jka=ikk-iqq
      jkk=mod(jka,ngkpt)
      do ii=1,3
        if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii)
      enddo
      ikptq=ikndx(jkk(1),jkk(2),jkk(3))
      isymq=isymndx(jkk(1),jkk(2),jkk(3))
      igsymq=isymg(:,ikptq,isymq)
      jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ictr
      igg=nint(dble(jkk-jka)/dble(ngkpt))
!write(6,*) '>>>>> ixt = ',ixt
!write(6,*) '>>>>> ixx = ',ixx
!write(6,*) '>>>>> iqq = ',iqq
!write(6,*) '>>>>> ikk = ',ikk
!write(6,*) '>>>>> jkk = ',jkk
!write(6,*) '>>>>> jka = ',jka
!write(6,*) '>>>>> igg = ',igg
!write(6,*) '>>>>> ikptq = ',ikptq,nkpt
!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt)
!write(6,'(a,3f10.5)') 'kptq = ',kpt(:,ikptq)
!read(*,*)
      qq=qq+ipwx(:,ipw1)
      kmq = kv - qq
      qq2=0.d0
      kmq2 = 0.d0
      do ii=1,3
      do jj=1,3
        qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj)
        kmq2=kmq2+kmq(ii)*bmet(ii,jj)*kmq(jj)
      enddo
      enddo
      vq(ix,iy,iz)=4.d0*pi/qq2
      ekmq = kmq2/(2.d0*pwse_mass)
!write(6,'(a,3f10.5)') 'qq = ',qq
!write(6,'(a,3f10.5)') 'qq2 = ',qq2
!write(6,*) 'vq = ',vq(ix,iy,iz)
!write(6,*) 'jka = ',jka
!write(6,*) 'jkk = ',jkk
!write(6,*) 'jkb = ',jkb
!write(6,*) 'igg = ',igg
!write(6,*) 'ikpt = ',ikpt
!write(6,*) 'ikptq = ',ikptq
      omega(ix,iy,iz)=ekmq-ek+dw/2
!write(33,'(3i3,2x,4f10.6)') ixt,ek,ekmq,ek-ekmq,9.879021/27.2114
      lqsing=.false.
      xmat2(ix,iy,iz)=(1.d0,0.d0)
      vqmat2(ix,iy,iz)=xmat2(ix,iy,iz)*vq(ix,iy,iz)
      if (mode.eq.1) vqmat2(ix,iy,iz) = vqmat2(ix,iy,iz)*(-omega(ix,iy,iz))
    enddo
    enddo
    enddo
!read(*,*)

! DEBUG
!write(6,*) "Matrix elements computed"

    if (itetrahedron.eq.0) then
      write(6,*) "Sum over points not yet implemented in self energy calculation"
      write(6,*) "set itetrahedron=1 and run again"
    else
! tetrahedron integration

! Don't waste time integrating over energies with no contribution
      do ix=1,nxtgkpt(1)
      do iy=1,nxtgkpt(2)
      do iz=1,nxtgkpt(3)
        ixt=(/ix,iy,iz/)
        ixx=ixt*ngkpt/nxtgkpt
        call fhilo(omega,ixt,nxtgkpt,whi,wlo)
        if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
          wwhi=whi
          wwlo=wlo
        else
          wwhi=max(wwhi,whi)
          wwlo=min(wwlo,wlo)
        endif
      enddo
      enddo
      enddo

      iwlo=nint(wwlo*dble(nwpt)/wmax)
      iwhi=nint(wwhi*dble(nwpt)/wmax)
      if (ocsign.lt.0.d0) then
        ielo=iwlo+1
        iehi=iwhi+nwpt
      else
        ielo=iwlo-nwpt
        iehi=iwhi-1
      endif
      allocate (ssi(ielo:iehi))
      do ie=ielo,iehi
        ssi(ie)=0.d0
      enddo
      ssx=0.d0

      if (lc) then
        do iw=1,nwpt
          do ix=1,nxtgkpt(1)
          do iy=1,nxtgkpt(2)
          do iz=1,nxtgkpt(3)
            ixt=(/ix,iy,iz/)
            ixx=ixt*ngkpt/nxtgkpt
            ixtqpt = ixtqndx(ixt(1),ixt(2),ixt(3))
            if (ipw1.le.npwx) then
              vfactor(ixt(1),ixt(2),ixt(3),iw)= &
&                 vqmat2(ixt(1),ixt(2),ixt(3))* &
&                 (xtlossfn(iw,ixtqpt,ipw1)-conjg(xtlossfn(iw,ixtqpt,ipw1)))*0.5d0
            else
              vfactor(ixt(1),ixt(2),ixt(3),iw)=(0.d0,0.d0)
            endif
!if (ipw1.eq.59.and.iw.eq.247.and.(ix.eq.7.or.ix.eq.8).and.(iy.eq.9.or.iy.eq.10).and.(iz.eq.7.or.iz.eq.8)) then
!write(34,*) ixt,vqmat2(ix,iy,iz),xtlossfn(iw,ixtqpt,jjpw)-conjg(xtlossfn(iw,ixtqpt,jjpw)),vfactor(ixt(1),ixt(2),ixt(3),iw)
!endif
!if (iw.eq.idum) write(33,'(3i3,2x,2f10.6,2x,f10.6)') ixt,xtlossfn(iw,ixtqpt,jjpw),omega(ix,iy,iz)
          enddo
          enddo
          enddo
!write(6,*) "iw = ", iw
        enddo
      endif

! DEBUG
!write(6,*) "Beginning integration"

! Do the integration
      do ix=1,nxtgkpt(1)
      do iy=1,nxtgkpt(2)
      do iz=1,nxtgkpt(3)
!xdum = 0.d0;
!xdum2 = 0.d0;
!xdum3 = 0.d0;
        ixt=(/ix,iy,iz/)
        ixx=ixt*ngkpt/nxtgkpt
        iqq=ixt-ixtctr
        ssx1=ssx
        call fval(omega,ixt,ixtp,nxtgkpt,enval)
        if (lqsing) call fval(vq,ixt,ixtp,nxtgkpt,vq2)
        if (lc) then
          do iw=1,nwpt
            call fpol(vfactor(:,:,:,iw),nxtgkpt,ixt,ixtp,vfv(:,iw))
          enddo
        endif
        if (lx) call fpol(vqmat2(:,:,:),nxtgkpt,ixt,ixtp,vfx)
        do itet=1,6
! DEBUG
!cdum = 0.d0
          do iv=1,4
            evtx0(iv)=enval(ivndx(iv,itet))
            rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet))
            if (lc) then
              do iw=1,nwpt
                vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw)
              enddo
            endif
            if (lx) xpyr0(iv)=vfx(ivndx(iv,itet))
          enddo
          call indxhpsort(4,4,evtx0,indxe)
          evtx=evtx0(indxe)
          do ii=1,4
            jj=indxe(ii)
            kvtx(:,ii) = (iqq(:) + rrpyr(:,jj))
            do kk=1,3
              kvtxc(kk,ii)=dot_product(iqq+rrpyr(:,jj),qkcvt(:,kk))
            enddo
          enddo
          xk(1:3,1)=(/0.d0,0.d0,0.d0/)
          xkc(1:3,1)=(/0.d0,0.d0,0.d0/)
          do ii=1,3
            xk(:,ii+1)=rrpyr(:,indxe(ii+1))-rrpyr(:,indxe(1))
          enddo
          do ii=2,4
            xkc(1:3,ii)=kvtxc(1:3,ii)-kvtxc(1:3,1)
          enddo
          do jj=1,3  ! make contragredient, reduced coordinates
            call cross(xk(:,mod(jj,3)+2),xk(:,mod(jj+1,3)+2),rg(:,jj))
          enddo
          do jj=1,3  ! make contragredient, Cartesian coordinates
            call cross(xkc(:,mod(jj,3)+2),xkc(:,mod(jj+1,3)+2),rgc(:,jj))
          enddo
          tvol=dot_product(xk(1:3,2),rg(1:3,1))
          rg=rg/tvol
          tvol=abs(tvol)
          tvolc=dot_product(xkc(1:3,2),rgc(1:3,1))
          rgc=rgc/tvolc
          tvolc=abs(tvolc)
          iwwlim=nint(evtx/dw)
          if (lc) then
            do iw=1,nwpt
              aa0(iw)=vpyr0(indxe(1),iw)   ! base value of function
              av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function
              do jj=1,3
                av(1:3,iw)=av(1:3,iw) &
&                      +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj)
              enddo
            enddo
            do iww=iwwlim(1)-1,iwwlim(4)+1
              wwwlo=dw*(iww-0.5d0)
!              wwwhi=dw*(iww+0.5d0)
              call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,wwwlo,dw,sint1,svec)
              do iw=1,nwpt
                ie=iww-ocsign*iw
                if (ie.ge.ielo.and.ie.le.iehi) ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt
              enddo
            enddo
          endif
! DEBUG
!write(32,'(4i3,2x,f14.10,2x,f14.10)') ixt,itet,(cdum/2.d0)
          if (lx) then
            xme=xpyr0(indxe(1))
            xv=(/0.d0,0.d0,0.d0/)
            do jj=1,3
              xv=xv+(xpyr0(indxe(jj+1))-xpyr0(indxe(1)))*rg(:,jj)
            enddo
            call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,evtx(1),evtx(4)-evtx(1),sint1,svec)
            ssx=ssx+(xme*sint1+dot_product(svec,xv))/volelmnt
          endif
!write(33,'(3i3,2x,i3,2x,4f10.6,2x,2f10.6)') ixt,itet,evtx,cdum
        enddo
      enddo
      enddo
      enddo
!    if (lc) ssi=-ssi/(dble(nxtgkpt(1)*nxtgkpt(2)*nxtgkpt(3))*vol*pi)
!    if (lx) ssx=-ssx/(dble(nxtgkpt(1)*nxtgkpt(2)*nxtgkpt(3))*vol)
      if (lc) ssi=-ocsign*ssi/((2*pi)**3*pi)
      if (lx) ssx=-ssx/((2*pi)**3)

! DEBUG
!write(6,*) "Integration complete"

      if (lc) then
        do ie=-nwpt,nwpt
          eps1=dble(ie)*dw-dw/2
!        ssc(ie)=(0.d0,0.d0)
          do je=ielo,iehi
            if (ie.eq.je) then
              ssc(ie)=ssc(ie)+pi*ssi(je)
            else
              eps2=dble(je)*dw-dw/2
              ssc(ie)=ssc(ie)-(0.d0,1.d0)*ocsign*ssi(je)*dw/(eps2-eps1)
            endif
!if (ie.eq.0) write(16,'(2i4,2f10.3,es12.3,2(2x,2es11.3))') ie,je,eps1*27.2114,eps2*27.2114,1/(eps2-eps1),ssc(ie)*27.2114,dble(ssi(je)*27.2114)
!if (ie.eq.0) write(16,'(2x,i4,2x,f10.3,4x,f10.6,2x,f10.6)') je,eps2*27.2114,ssi(je)*27.2114
          enddo
        enddo
      endif
!do je=ielo,iehi
!write(6,*) je,ssi(je)
!enddo
      deallocate (ssi)
!write(6,'("lstv",2x,f14.10,4x,"(",f14.10,",",f14.10,")")') ssx,(ssc(0)+ssc(1))/2.

      if (lc) cse1=cse1+ssc
      if (lx) xse=xse+ssx
    endif

  enddo

if (itetrahedron.eq.1) then
  cse=(cse1(0)+cse1(1))/2
  dcse=(cse1(1)-cse1(0))/dw
endif
zz=1.d0/(1.d0-dcse)
!write(6,*) cse*27.2114
!write(6,*) xse*27.2114
!write(6,*) zz
!write(6,*) dble(zz)*dimag(cse)*27.2114

!temparray=log10(temparray)
!temparray=10**(temparray-int(temparray))
!do ix=1,ngkpt(1)
!  write(76,*) ix,'--------------------------------------------------'
!  do iy=1,ngkpt(2)
!!    write(76,'(10es8.1)') dble(temparray(ix,iy,:))
!    write(76,'(10f8.5)') dble(temparray(ix,iy,:))
!  enddo
!enddo

!write(6,*) "AAA2", time()-xdum

return
end subroutine mkxtpwse

!*************************************************************************

double precision function re_dielf_feg(omega)
implicit none
double precision :: omega
double precision :: qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi
double precision :: NN
double precision :: EE,GG,hh,pf,AA,BB,Elo,Ehi
double precision :: lnterm1,lnterm2,lnterm3,lnterm4,lnterm5
double precision :: xterm1,xterm2,xterm3,xterm4,xterm5,rterm
common /feg/ qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi

  EE = omega/Ef
  GG = Eg/Ef
  hh = qq/xkf
  pf = (3.d0/4.d0)*(omegap2/(Ef*Ef))
  AA = GG + hh*hh
  BB = 0.25/(hh*hh)
  Ehi = GG + hh*abs(hh+2)
  Ehi = GG + hh*abs(hh-2)
  if (hh.le.2.d0) then
    NN = (3*GG+4*hh-0.25*GG*hh*hh)*hh*hh
    rterm = 2*(Elo-GG) - 4*hh*hh*(GG+2*hh-2*AA)*BB
  else
    NN = 4*hh*(GG+hh*hh)
    rterm = 8*hh*AA*BB
  endif
  xterm1 = 1.d0 - (AA+EE)**2/BB
  xterm2 = 1.d0 - (AA-EE)**2/BB
  lnterm1 = log(abs((Ehi+EE)/(Elo+EE)))
  lnterm2 = log(abs((Ehi-EE)/(Elo-EE)))
  re_dielf_feg = 1.d0 + (pf/NN)*(rterm + xterm1*lnterm1 + xterm2*lnterm2)
  if (hh.le.2.d0) then
    xterm3 = -(GG+EE)
    xterm4 = -(GG-EE)
    lnterm3 = log(abs((Elo+EE)/(GG+EE)))
    lnterm4 = log(abs((Elo-EE)/(GG-EE)))
    re_dielf_feg = re_dielf_feg + (pf/NN)*(xterm3*lnterm3 + xterm4*lnterm4)
  endif

end function re_dielf_feg

!*************************************************************************

double precision function im_dielf_feg(omega)
implicit none
double precision :: omega
double precision :: qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi
double precision :: term,NN
double precision :: EE,GG,hh,pf
common /feg/ qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi

!write(*,*) edens,Eg,omegap2
!read(*,*)
  EE = omega/Ef
  GG = Eg/Ef
  hh = qq/xkf
  pf = (3.d0*pi/4.d0)*(omegap2/(Ef*Ef))
!write(*,*) pf,omegap2,Ef
!read(*,*)
  if (hh.le.2.d0) then
    NN = (3*GG+4*hh-0.25*GG*hh*hh)*hh*hh
    if (EE.le.GG) then
      im_dielf_feg = 0.d0
    else if (EE.le.(GG+hh*(2.d0-hh))) then
      im_dielf_feg = pf * (EE-GG)/NN
    else if (EE.le.(GG+hh*(2.d0+hh))) then
      term = (EE-GG-hh*hh)/(2*hh)
      im_dielf_feg = pf * (1.d0 - term**2)/NN
    else
      im_dielf_feg = 0.d0
    endif
  else
    NN = 4*hh*(GG+hh*hh)
    if (EE.le.(GG+hh*(hh-2.d0))) then
      im_dielf_feg = 0.d0
    else if (EE.le.(GG+hh*(2.d0+hh))) then
      term = (EE-GG-hh*hh)/(2*hh)
      im_dielf_feg = pf * (1.d0 - term**2)/NN
    else
      im_dielf_feg = 0.d0
    endif
  endif

end function im_dielf_feg

!*************************************************************************

double complex function dielf_feg(omega)
implicit none
double precision :: omega
double precision :: qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi
double precision :: eps1, eps2
double precision :: re_dielf_feg,im_dielf_feg
external re_dielf_feg,im_dielf_feg
common /feg/ qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi

  eps1 = re_dielf_feg(omega)
  eps2 = im_dielf_feg(omega)+broadening
  dielf_feg = eps1 + (0.d0,1.d0)*eps2

end function dielf_feg

!*************************************************************************

double precision function dre_dielf_feg(omega)
implicit none
double precision :: omega
double precision :: qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi
double precision :: lnterm1,lnterm2,xterm1,xterm2,xterm3,xterm4
double precision :: dlnterm1,dlnterm2
common /feg/ qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi

  lnterm1 = log((Eq+qkf+omega)/(Eq-qkf+omega))
  lnterm2 = log((Eq+qkf-omega)/(Eq+qkf+omega))
  dlnterm1 = 1.d0/(Eq+qkf+omega)-1.d0/(Eq-qkf+omega)
  dlnterm2 = 1.d0/(Eq+qkf-omega)-1.d0/(Eq+qkf+omega)
  xterm1 = (4*Ef*Eq-(Eq+omega)**2)*dlnterm1/(2*qkf*(qq**2))
  xterm2 = (4*Ef*Eq-(Eq-omega)**2)*dlnterm2/(2*qkf*(qq**2))
  xterm3 = -(Eq+omega)*lnterm1/(2*qkf*(qq**2))
  xterm4 = (omega-Eq)*lnterm2/(2*qkf*(qq**2))
  dre_dielf_feg = 0.5d0*(qtf/qq)**2*(xterm1+xterm2+xterm3+xterm4)

end function dre_dielf_feg

!*************************************************************************

subroutine lossfn_free_electron_gas(omegahi,omegalo,omegap,qq0,broadening0,Ef0,qtf0,xkf0,edens0,lossfn_feg)
implicit none
double precision :: omegahi,omegalo,omegap
double precision :: qq0,Ef0,qtf0,xkf0,edens0,broadening0
double precision :: qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi
double precision :: abr,rlr,error
integer :: nsing,numcal,maxns
double precision :: xsing(20),omegaq,wtest
double complex :: lossfn_feg
double complex :: dielf_feg,cgrater
external dielf_feg,cgrater
common /feg/ qq,Eq,Ef,qtf,xkf,qkf,edens,broadening,Eg,omegap2,pi

  qq = qq0
  Ef = Ef0
  qtf = qtf0
  xkf = xkf0
  edens = edens0
  Eq = 0.5d0*qq*qq
  qkf = qq*xkf
  broadening = broadening0
  omegaq = omegap + 0.3*qkf*qkf/omegap

  abr = 1.d-7
  rlr = 1.d-7

  nsing = 0
  wtest = omegaq-100.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq-10.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq-3.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq-1.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq+1.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq+3.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq+10.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif
  wtest = omegaq+100.*broadening
  if ((omegahi.gt.wtest).and.(wtest.gt.omegalo)) then
    nsing = nsing+1
    xsing(nsing) = wtest
  endif

  lossfn_feg = cgrater(dielf_feg,omegalo,omegahi,abr,rlr,nsing,xsing,error,numcal,maxns)/(omegahi-omegalo)

end subroutine lossfn_free_electron_gas

max_stars_repo_name: miroi/nwchem

max_issues_repo_name: miroi/nwchem

max_forks_repo_name: miroi/nwchem


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! test host's memory arguments in an internal procedure
!

module foo
  contains
    subroutine sub(result,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14)
      integer result(:)
      integer d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14
      integer l1
      l1 = 25
      call csub
    contains
      subroutine csub()
	result(1) = d1
	result(2) = d2
	result(3) = d3
	result(4) = d4
	result(5) = d5
	result(6) = d6
	result(7) = d7
	result(8) = d8
	result(9) = d9
	result(10) = d10
	result(11) = d11
	result(12) = d12
	result(13) = d13
	result(14) = d14
	result(15) = l1
      end subroutine
    end subroutine
end module

program p
use foo
  integer result(15)
  integer expect(15)
  data expect/ 11,12,13,14,15,16,17,18,19,20,21,22,23,24,25/

  call sub(result,11,12,13,14,15,16,17,18,19,20,21,22,23,24)
  call check(result,expect,15)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
Kober Apartments is...

To learn more about rental housing in Davis, check out our Housing Guide as well as our Apartments pages.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine r8mkspl3pb(fun,x,nx,th,nth,ph,nph,fspl,nf4,nf5,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,nb1,
     >   wk,inwk,ilinx,ilinth,ilinph,ier)
C
C  create a tricubic biperiodic spline with knots at grid points and
C  function values from the callable function `fun' passed.
C
C  use tpsplinb to setup the spline coefficients
C
C  periodic boundary condition for theta & phi grids.
C  boundary condition data may be specified at x(1) and x(nx) for each
C  theta & phi:
C     ibcxmin=0 ==> "not a knot" boundary condition, cf cubspl.for
C     ibcxmin=1 ==> match slope, bcxmin(ith,iph) gives df/dx at
C         x(1),th(ith),ph(iph)
C     ibcxmin=2 ==> match 2nd derivative, given at x(1),th(ith),ph(iph)
C          by bcxmin(ith,iph)
C
C     ibcxmax,bcxmax have analogous interpretation -- at x(nx).
C
!============
! idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100)
      INTEGER nth,nph,nf4,nf5,nb1,inwk,nx,iph,ith,ix
!============
! idecl:  explicitize implicit REAL declarations:
      REAL*8 fun
!============
      external fun                      ! passed real function(x,th)
      REAL*8 x(nx)                        ! x coordinate array
      REAL*8 th(nth)                      ! th coordinate array
      REAL*8 ph(nph)                      ! ph coordinate array
C
      REAL*8 fspl(4,4,4,nf4,nf5,nph)      ! function data / spline coeff array
      REAL*8 wk(inwk)                     ! workspace for bpsplinb
C
      integer ibcxmin                   ! boundary condition indicator
      REAL*8 bcxmin(nb1,nph)              ! boundary condition data
      integer ibcxmax                   ! boundary condition indicator
      REAL*8 bcxmax(nb1,nph)              ! boundary condition data
C
      integer ilinx                     ! output =1 if x(...) evenly spaced
      integer ilinth                    ! output =1 if th(...) evenly spaced
      integer ilinph                    ! output =1 if ph(...) evenly spaced
C
      integer ier                       ! completion code from bpspline
C
C----------------------------
C
      ier=0
      if(nf4.lt.nx) then
         write(6,'('' ?mkspl3pb -- array dim error, nf4 .lt. nx'')')
         ier=1
      endif
      if(nf5.lt.nth) then
         write(6,'('' ?mkspl3pb -- array dim error, nf5 .lt. nth'')')
         ier=2
      endif
      if(nb1.lt.nth) then
         write(6,'('' ?mkspl3pb -- array dim error, nb1 .lt. nth'')')
         ier=3
      endif
      if(inwk.lt.5*max(nx,nth,nph)) then
         write(6,'('' ?mkspl3pb -- array dim error, inwk too small'')')
         ier=4
      endif

      if(ier.ne.0) return
C
      do iph=1,nph
         do ith=1,nth
            do ix=1,nx
               fspl(1,1,1,ix,ith,iph)=fun(x(ix),th(ith),ph(iph))
            enddo
         enddo
      enddo
C
      call r8tpsplinb(x,nx,th,nth,ph,nph,fspl,nf4,nf5,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,nb1,
     >   wk,inwk,ilinx,ilinth,ilinph,ier)
C
      return
      end

max_stars_repo_name: ORNL-Fusion/RFSciDAC-testing

max_issues_repo_name: ORNL-Fusion/RFSciDAC-testing

max_forks_repo_name: ORNL-Fusion/RFSciDAC-testing


module DataInputM
  use UtilitiesM
  use DebuggerM

  use NodeM
  use NodePtrM
  use GeometryM
  use SourceM
  use ThermalMaterialM
  use ThermalElementM
  use Thermal2DApplicationM
  use ConvectionOnLineM
  use FluxOnLineM

  use MeshM
  use ModelM

  implicit none

  private
  public :: initFEM2D

  integer(ikind), parameter    :: projectData = 1
  integer(ikind), parameter    :: project = 2
  integer(ikind), dimension(8) :: date_time
  integer(ikind)               :: nElem
  integer(ikind)               :: nTriangElem
  integer(ikind)               :: nRectElem
  integer(ikind)               :: nPoint
  integer(ikind)               :: iPoint
  integer(ikind)               :: nNormalFlux
  integer(ikind)               :: nDirichlet
  integer(ikind)               :: nMaterial
  integer(ikind)               :: nGauss
  integer(ikind)               :: nConvection
  integer(ikind)               :: isQuadratic
  integer(ikind)               :: nSourceOnPoints
  integer(ikind)               :: nSourceOnSurfaces
  integer(ikind)               :: nPointSource
  integer(ikind)               :: nSurfaceSource
  integer(ikind)               :: printStep
  real(rkind)                  :: t0
  real(rkind)                  :: errorTol
  character(100)               :: projectName
  character(100)               :: aux
  logical       , parameter    :: verbose = .false.

  interface initFEM2D
     procedure :: initFEM2D
  end interface initFEM2D

contains

  subroutine initFEM2D(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    print'(A)', 'Initializing Thermal2D application'
    call initLog(.true., 'log.dat')
    call debugLog('  Reading project data')
    call readProjectData
    call debugLog('  Reading mesh data')
    call initMesh(thermalAppl)
    call debugLog('  Reading materials properties')
    call initMaterials(thermalAppl)
    call debugLog('  Reading elements')
    call initElements(thermalAppl)
    call debugLog('  Reading point and line Sources')
    call readPointLineSurfaceSources(thermalAppl)
    call debugLog('  Reading Boundary Conditions')
    call readBoundaryConditions(thermalAppl)
    !call debugLog('  Reading Initial Values')
    !call readInitialValues(thermalAppl)
    call debugLog('End loading data')
  end subroutine initFEM2D

  subroutine readProjectData
    implicit none
    open(projectData, file = 'projectData.dat')
    read(projectData, '(*(A))') projectName
    close(projectData)
    call debugLog('    Project name: ', trim(projectName))
  end subroutine readProjectData

  subroutine initMesh(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    integer(ikind) :: i
    real(rkind)    :: x, y, z
    open(project, file = trim(projectName)//'.dat')
    do i = 1, 9
       read(project,*)
    end do
    read(project,*)  aux, nElem
    read(project,*)  aux, nPoint
    read(project,*)  aux, isQuadratic
    read(project,*)  aux, nTriangElem
    read(project,*)  aux, nRectElem
    read(project,*)  aux, nMaterial
    read(project,*)  aux, nGauss
    read(project,*)  aux, nDirichlet
    read(project,*)  aux, nNormalFlux
    read(project,*)  aux, nConvection
    read(project,*)  aux, nSourceOnPoints
    read(project,*)  aux, nSourceOnSurfaces
    read(project,*)  aux, nPointSource
    read(project,*)  aux, nSurfaceSource
    read(project,*)  aux, printStep
    read(project,*)  aux, t0
    read(project,*)  aux, errorTol
    call debugLog('    Number of Elements.............................: ', nElem)
    call debugLog('    Are Elements Quadratic.........................: ', isQuadratic)
    call debugLog('    Number of Triangular elements..................: ', nTriangElem)
    call debugLog('    Number of Rectangular elements.................: ', nRectElem)
    call debugLog('    Number of Nodes................................: ', nPoint)
    call debugLog('    Number of Dirichlet conditions.................: ', nDirichlet)
    call debugLog('    Number of NormalFluxOnLines conditions.........: ', nNormalFlux)
    call debugLog('    Number of ConvectionOnLines conditions.........: ', nConvection)
    call debugLog('    Number of Sources on points....................: ', nSourceOnPoints)
    call debugLog('    Number of Sources on surfaces..................: ', nSourceOnSurfaces)
    call debugLog('    Number of points with pointSource..............: ', nPointSource)
    call debugLog('    Number of Surfaces with surfaceSource..........: ', nSurfaceSource)
    call debugLog('    Number of Materials............................: ', nMaterial)
    call debugLog('    Gauss cuadrature order.........................: ', nGauss)
    call debugLog('    PrintStep......................................: ', printStep)
    call debugLog('    Initial time...................................: ', t0)
    call debugLog('    Error tolerance................................: ', errorTol)

    thermalAppl = thermal2DApplication(                        &
           nNode = nPoint                                      &
         , nElement = nTriangElem + nRectElem                  &
         , nNormalFlux = nNormalFlux                           &
         , nConvection = nConvection                           &
         , nSource = nSourceOnPoints + nSourceOnSurfaces       &
         , nMaterial = nMaterial                               &
         , nGauss = nGauss                                     )

    call thermalAppl%setTransientValues(printStep, t0, errorTol)

    do i = 1, 6
       read(project,*)
    end do
    if(verbose) print'(A)', 'Coordinates:'
    if(verbose) print'(1X,A)', 'NODES      X          Y'
    do i = 1, nPoint
       read(project,*) iPoint, x, y
       if(verbose) print'(3X,I0,3X,E10.3,3X,E10.3,3X,E10.3)', iPoint, x, y
       thermalAppl%node(iPoint) = node(iPoint, 1, x, y)
       call thermalAppl%node(iPoint)%assignDof(1, thermalAppl%model%dof(iPoint))
       call thermalAppl%model%addNode(iPoint, thermalAppl%node(iPoint))
    end do
  end subroutine initMesh

  subroutine initMaterials(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    integer(ikind) :: i, iMat
    real(rkind) :: kx, ky, cp, rho
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(A)', 'Material         Kx            Ky           cp          rho       '
    do i = 1, nMaterial
       read(project,*) iMat, kx, ky, cp, rho
       thermalAppl%material(iMat) = thermalMaterial(kx, ky, cp, rho)
       if(verbose) print'(4X,I0,7X,2(E10.3,3X))', iMat, kx, ky, cp, rho
    end do
  end subroutine initMaterials

  subroutine initElements(thermalAppl)
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    type(NodePtrDT), dimension(:), allocatable :: auxNode
    integer(ikind) :: i, j, iElem, iMat, nNode, Conectivities(8)
    character(len=13) :: type
    Conectivities = 0
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(A)', 'Element  |      Type      |  material index  |  nNodes  |  connectivities'
    do i = 1, nElem
       read(project,*) iElem, type, iMat, nNode, (Conectivities(j),j=1,nNode)
       if(verbose) print'(I5,A15,I18,I14,5X,*(I5,X))', iElem, type, iMat, nNode, (Conectivities(j),j=1,nNode)
       allocate(auxNode(nNode))
       do j = 1, nNode
          call auxNode(j)%associate(thermalAppl%node(conectivities(j)))
       end do
       thermalAppl%element(iElem) = thermalElement(iElem, auxNode, thermalAppl%material(iMat))
       call thermalAppl%model%addElement(i, thermalAppl%element(iElem))
       deallocate(auxNode)
    end do
  end subroutine initElements

  subroutine readPointLineSurfaceSources(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    integer(ikind)                              :: i, countSource, auxInt
    integer(ikind)                              :: iNode, iElem, iSource
    character(150), dimension(1)                :: func
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'nSource'
    if(verbose) print'(A)', 'Source    Function'
    do i = 1, nSourceOnPoints+nSourceOnSurfaces
       read(project,*) iSource, func(1)
       thermalAppl%source(iSource) = source(2, 1, (/'x', 'y'/), func)
       if(verbose) print'(I0,5X,A)', iSource, func(1)
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'pointSources'
    if(verbose) print'(A)', 'Node    Source'
    do i = 1, nPointSource
       read(project,*) iNode, iSource
       if(verbose) print'(I0,5X,I0)', iNode, iSource
       call thermalAppl%node(iNode)%assignSource(thermalAppl%source(iSource))
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'surfaceSources'
    if(verbose) print'(A)', 'Element   Source'
    do i = 1, nSurfaceSource
       read(project,*) iElem, iSource
       if(verbose) print'(I0,5X,I0)', iElem, iSource
       call thermalAppl%element(iElem)%assignSource(thermalAppl%source(iSource))
    end do
  end subroutine readPointLineSurfaceSources

  subroutine readBoundaryConditions(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout)  :: thermalAppl
    integer(ikind)                               :: i, j, id, elemID, nPointID
    integer(ikind)                               :: iPoint, conditionCounter
    integer(ikind), dimension(:), allocatable    :: pointID
    real(rkind)                                  :: value
    real(rkind)                                  :: coef, temp
    type(ThermalElementDT)                       :: element
    type(NodePtrDT)  , dimension(:), allocatable :: node
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Dirichlet conditions'
    if(verbose) print'(A)', 'Node    Value'
    do i = 1, nDirichlet
       read(Project,*) id, value
       if(verbose) print'(I0,5X,E10.3)', id, value
       call thermalAppl%node(id)%fixDof(1, value)
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(isQuadratic == 0) then
       nPointID = 2
    else
       nPointID = 3
    end if
    allocate(pointID(nPointID))
    allocate(node(nPointID))
    conditionCounter = 0
    if(verbose) print'(/,A)', 'Normal Flux On Lines conditions'
    if(verbose) print'(A)', 'Elem    Nodes     Value'
    do i = 1, nNormalFlux
       conditionCounter = conditionCounter + 1
       read(Project,*) elemID, (pointID(j),j=1,nPointID), value
       if(verbose) print*, elemID, (pointID(j),j=1,nPointID), value
       element = thermalAppl%element(elemID)
       do j = 1, nPointID
          node(j) = element%node(pointID(j))
       end do
       thermalAppl%normalFluxOL(i) = fluxOnLine(i, pointID, value, node, element%geometry)
       call thermalAppl%model%addCondition(conditionCounter, thermalAppl%normalFluxOL(i))
    end do
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Convection On Lines conditions'
    if(verbose) print'(A)', 'Elem    Nodes     Coef     Temp'
    do i = 1, nConvection
       conditionCounter = conditionCounter + 1
       read(Project,*) elemID, (pointID(j),j=1,nPointID), coef, temp
       if(verbose) print*, elemID, (pointID(j),j=1,nPointID), coef, temp
       element = thermalAppl%element(elemID)
       do j = 1, nPointID
          node(j) = element%node(pointID(j))
       end do
       thermalAppl%convectionOL(i) = convectionOnLine(i, pointID, coef, temp, node, element%geometry)
       call thermalAppl%model%addCondition(conditionCounter, thermalAppl%convectionOL(i))
    end do
  end subroutine readBoundaryConditions

  subroutine readInitialValues(thermalAppl)
    implicit none
    type(Thermal2DApplicationDT), intent(inout) :: thermalAppl
    integer(ikind) :: i, nodeID
    real(rkind) :: temp
    do i = 1, 7
       read(project,*)
    end do
    if(verbose) print'(/,A)', 'Initial Temperature'
    if(verbose) print'(A)', 'Node   Temperature'
    do i = 1, nPoint
       read(project,*) nodeID, temp
       if(verbose) print*, nodeID, temp
       thermalAppl%model%dof(nodeID) = temp
    end do
  end subroutine readInitialValues

end module DataInputM


      SUBROUTINE dcosft2d(Y,n1_max,N1,n2,wrk)
!
! This is a fast cosine transform routine for a 2-d function Y(n1,n2).
! n1-1 and n2-1 must both be powers of 2 and >= 2.
! This function is it's own transform, but the inverse must be
! multiplied by 4/((n1-1)*(n2-1)).   Since the data and transform are
! all real, we simply FCT in x along the rows, then in y along the
! columns using the 1-D routine.
!
      IMPLICIT NONE
      INTEGER n1_max, n1, n2
      REAL*8  Y(N1_max,n2), wrk(n2)
      INTEGER j, k

      DO k=1,n2
         CALL dcosft(y(1,k),n1-1)
      END DO
      DO j=1,n1
         DO k=1,n2
            wrk(k)=y(j,k)
         END DO
         CALL dcosft(wrk,n2-1)
         DO k=1,n2
            y(j,k)=wrk(k)
         END DO
      END DO

      RETURN
      END SUBROUTINE DCOSFT2D


text:
      subroutine saxpy(sx,sy)

      INCLUDE 'SAXPY_INC.f'

c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c

      real sx(SIZE),sy(SIZE),sa
      integer i,incx,incy,ix,iy,m,mp1

      sa = 3
      incx = 1
      incy = 1
c
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      do 10 i = 1,SIZE
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c

c 20   PRINT *, 'bef'
 20   do 50 i = 1,SIZE
        sy(i) = sy(i) + sa*sx(i)
   50 continue
c      PRINT *, 'aft'
      return
      return
      end

      subroutine saxpyref(sx,sy)

      INCLUDE 'SAXPY_INC.f'

c
c     constant times a vector plus a vector.
c     uses unrolled loop for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
      real sx(SIZE),sy(SIZE),sa
      integer i,incx,incy,ix,iy,m,mp1

      sa = 3
      incx = 1
      incy = 1
c
      if (sa .eq. 0.0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      do 10 i = 1,SIZE
        sy(iy) = sy(iy) + sa*sx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c

 20   i=1
      do 50 i = 1,SIZE
        sy(i) = sy(i) + sa*sx(i)
   50 continue
c      PRINT *, 'aft2'
      return
      return
      end

      INTEGER FUNCTION CLOCK()

      CLOCK = 0

      RETURN
      END

      PROGRAM MAIN

      INCLUDE 'SAXPY_INC.f'

      INTEGER i, j, incx, incy
      INTEGER bef, af, diff
      INTEGER befRef, afRef, diffRef
      INTEGER CLOCK
      REAL  sa, sx(1:SIZE), sy(1:SIZE)
      REAL  sxRef(1:SIZE), syRef(1:SIZE)

      LOGICAL success

      incx = 1
      incy = 1
      sa = 3
      DO 10 i = 1, SIZE
          sx(i) = i
          sy(i) = i
          sxRef(i) = i
          syRef(i) = i
10    ENDDO

      bef = CLOCK()
      DO 20 i = 1, 10000
          CALL saxpy(sx,sy)
20    ENDDO
      af = CLOCK()

      befRef = CLOCK()
      DO 21 i = 1, 10000
          CALL saxpyRef(sxRef,syRef)
21    ENDDO
      afRef = CLOCK()

      success = .TRUE.

      DO 22 i = 1, SIZE
         IF((ABS(sy(i) - syRef(i)) / syRef(i)) .GT. 0.000001) THEN
             success = .FALSE.
         ENDIF
22    ENDDO
      PRINT *, sy(100), syRef(100)
      IF(success .EQV. .TRUE.) THEN
         PRINT *, 'SUCCESS'
      ELSE
         PRINT *, 'ERROR'
      ENDIF

      diff = af - bef
      diffRef = afRef - befRef
      PRINT *, 'time: ', diff
      PRINT *, 'reference time: ', diffRef

      RETURN
      END


text:
! @@name:	atomic_restrict.2f
! @@type:	F-fixed
! @@compilable:	maybe
! @@linkable:	no
! @@expect:	failure
      SUBROUTINE SUB()
        COMMON /BLK/ R
        REAL R

!$OMP   ATOMIC UPDATE
          R = R + 1.0
      END SUBROUTINE SUB

      SUBROUTINE ATOMIC_WRONG2()
        COMMON /BLK/ I
        INTEGER I

!$OMP   PARALLEL

!$OMP     ATOMIC UPDATE
            I = I + 1
          CALL SUB()
!$OMP   END PARALLEL
      END SUBROUTINE ATOMIC_WRONG2

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


! chk_submodule.f90
!     Check: does the compiler support sumodules?
!
module base_module
    implicit none

    interface
        real module function add_some_number( value )
            real, intent(in) :: value
        end function add_some_number
    end interface
end module base_module

submodule(base_module) actual_function
contains
module procedure add_some_number
    some_number = 1.0 + value
end procedure add_some_number
end submodule actual_function

program chk_submodule
    use base_module

    implicit none

    real :: x, y

    x = 1.0
    y = add_some_number( x )

    write( *, '(a)'   ) 'Submodules are supported'

end program chk_submodule