text:
MODULE xc_vdw_module

  use grid_module, only: grid
  use xc_variables, only: xcene, xcpot
  use ggrid_module, only: NMGL, MGL, GG, NGgrid, LLG, MG_0, MG_1 &
                         ,construct_ggrid, destruct_ggrid
  use gradient_module
  use parallel_module, only: comm_grid, myrank, nprocs, ir_grid, id_grid
  use fd_module, only: fd, construct_nabla_fd, destruct_nabla_fd
  use lattice_module, only: lattice,get_aa_lattice,get_reciprocal_lattice
  use basic_type_factory
  use fft_module
  use polint_module

  implicit none

  PRIVATE
  PUBLIC :: init_xc_vdw, calc_xc_vdw

  integer,parameter :: DP=kind(0.0d0)
  integer,parameter :: QP=kind(0.0q0)

  real(8) :: Qmax=5.0d0
  real(8) :: Qmin=0.0d0
  integer :: NumQgrid=20
  real(8) :: SizeQgrid
  real(8),allocatable :: Qgrid(:)
  real(8) :: Dmax=48.6d0
  real(8),parameter :: Zab=-0.8491d0
  integer :: Mmax=12

  real(8),allocatable :: phiG(:,:,:)
  real(8),allocatable :: bmat(:,:),cmat(:,:),dmat(:,:)

  integer :: NG_table=1000
  real(8) :: Gmax_table=50.0d0
  real(8),allocatable :: phiG_table(:,:,:)

  real(8),parameter :: zero_density = 1.d-10
  logical :: flag_init=.false.
  real(8) :: Ex_LDA, Ec_LDA, Ec_vdw

  integer,parameter :: unit_phig=9

  include 'mpif.h'

CONTAINS

  SUBROUTINE init_xc_vdw
    implicit none
    integer :: i,j,n,np,info
    logical :: file_exist

    if ( flag_init ) return

    call write_border( 0, " init_xc_vdw(start)" )

    SizeQgrid = Qmax / NumQgrid
    allocate( Qgrid(0:NumQgrid) ) ; Qgrid=0.0d0
    do i=1,NumQgrid
       Qgrid(i) = Qmin + i*SizeQgrid
    end do

    allocate( phiG(NMGL,0:NumQgrid,0:NumQgrid) )
    phiG=0.0d0

    allocate( phiG_table(NG_table,0:NumQgrid,0:NumQgrid) )
    phiG_table=0.0d0

   !call read_file_phiG( file_exist, phiG )
    call read_file_phiG( file_exist, phiG_table )

    if ( .not.file_exist ) then

       np=( NumQgrid*(NumQgrid+1) )/2
       n=0
       do j=1,NumQgrid
       do i=j,NumQgrid

          n = n + 1
          if ( np-1 < myrank .or. mod(n-1,nprocs) /= myrank ) cycle

         !call calc_phiG( Qgrid(i), Qgrid(j), phiG(:,i,j) )
         !if ( i /= j ) phiG(:,j,i) = phiG(:,i,j)
          call calc_phiG_table( Qgrid(i), Qgrid(j), phiG_table(:,i,j) )
          if ( i /= j ) phiG_table(:,j,i) = phiG_table(:,i,j)

       end do ! i
       end do ! j

      !call MPI_ALLREDUCE( MPI_IN_PLACE, phiG, size(phiG), MPI_REAL8 &
      !                   ,MPI_SUM, MPI_COMM_WORLD, info )
       call MPI_ALLREDUCE( MPI_IN_PLACE, phiG_table, size(phiG_table) &
                          ,MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, info )

      !call write_file_phiG( phiG )
       call write_file_phiG( phiG_table )

    end if

    call make_phiG_by_interpolation( phiG_table, phiG )

    deallocate( phiG_table )

    call init1_xc_vdw ! coefficients of 3rd spline

    flag_init = .true.

    call write_border( 0, " init_xc_vdw(end)")

  END SUBROUTINE init_xc_vdw

  SUBROUTINE read_file_phiG( flag, phiGread )
    implicit none
    logical,intent(OUT) :: flag
    real(8),intent(INOUT) :: phiGread(:,0:,0:)
    integer :: i,j,k,m,n
    flag=.false.
    if ( myrank == 0 ) then
       inquire( file="xc_vdw_phig.dat", exist=flag )
       if ( flag ) then
          open(unit_phig,file="xc_vdw_phig.dat",status="old")
          read(unit_phig,*) m,n
          if ( m == size(phiGread,1) .and. n == NumQgrid ) then
             write(*,*) "phiG is read from 'xc_vdw_phig.dat'"
             do k=1,n
             do j=k,n
                do i=1,m
                   read(unit_phig,*) phiGread(i,j,k)
                end do
                if ( j /= k ) phiGread(:,k,j)=phiGread(:,j,k)
             end do ! j
             end do ! k
          else
             flag=.false.
          end if
          close(unit_phig)
       end if
    end if
    call MPI_BCAST( flag, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, i )
    if ( flag ) then
       call MPI_ALLREDUCE( MPI_IN_PLACE, phiGread, size(phiGread) &
                          ,MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, i )
    end if
  END SUBROUTINE read_file_phiG

  SUBROUTINE write_file_phiG( phiGout )
    implicit none
    real(8),intent(IN) :: phiGout(:,0:,0:)
    integer :: i,j,k
    if ( myrank == 0 ) then
       write(*,*) "phiG is written to 'xc_vdw_phig.dat'"
       open(unit_phig,file="xc_vdw_phig.dat")
       write(unit_phig,*) size(phiGout,1), NumQgrid
       do k=1,NumQgrid
       do j=k,NumQgrid
          do i=1,size(phiGout,1)
             write(unit_phig,*) phiGout(i,j,k)
          end do
       end do ! j
       end do ! k
       close(unit_phig)
    end if
  END SUBROUTINE write_file_phiG

  SUBROUTINE calc_phiG_table( qi, qj, phiGtab )
    implicit none
    real(8),intent(IN)  :: qi,qj
    real(8),intent(OUT) :: phiGtab(:)
    integer :: nr,i
    real(8) :: rmax,dr,r,dG,G
    real(8),allocatable :: phi(:)

    phiGtab = 0.0d0

    rmax = Dmax/max( qi, qj )
    nr   = 400
    dr   = rmax/nr
    dG   = Gmax_table/NG_table

    allocate( phi(0:nr) ) ; phi=0.0d0

    do i=1,nr
       r = i*dr
       call calc_phi_b( qi*r, qj*r, phi(i) )
    end do ! i

    do i=1,NG_table
       G=(i-1)*dG
       call calc_phiG_sub( nr, dr, phi(1), G, phiGtab(i) )
    end do

    deallocate( phi )

  END SUBROUTINE calc_phiG_table

  SUBROUTINE calc_phiG_sub( nr, dr, phi, G, phiG )
    implicit none
    integer,intent(IN) :: nr
    real(8),intent(IN) :: dr, phi(nr), G
    real(8),intent(OUT) :: phiG
    integer :: i,ig
    real(8) :: r,sum0,pi4

    pi4 = 4.0d0*acos(-1.0d0)

    if ( G <= 1.d-9 ) then

       sum0=0.0d0
       do i=1,nr
          r = i*dr
          sum0 = sum0 + r*r*phi(i)
       end do
       phiG=pi4*sum0*dr

    else

       sum0=0.0d0
       do i=1,nr
          r = i*dr
          sum0 = sum0 + r*phi(i)*sin(G*r)
       end do
       phiG=pi4/G*sum0*dr

    end if

  END SUBROUTINE calc_phiG_sub

  SUBROUTINE calc_phiG( qi, qj, phiG )
    implicit none
    real(8),intent(IN)  :: qi,qj
    real(8),intent(OUT) :: phiG(NMGL)
    integer :: nr,ig,i
    real(8) :: rmax,dr,r,pi4,G,sum0,ct0,ct1,sum1
    real(8),allocatable :: phi(:), phiG0(:)

    phiG = 0.0d0

    rmax = Dmax/max( qi, qj )
    nr   = 400
    dr   = rmax/nr

    allocate( phi(0:nr)   ) ; phi=0.0d0
    allocate( phiG0(NMGL) ) ; phiG0=0.0d0

    do i=1,nr

       r = i*dr
       call calc_phi_b( qi*r, qj*r, phi(i) )

       call check_phiG( i, dr, phi(1), phiG )
       sum0=sum((phiG-phiG0)**2)
       if ( sum0 <= 1.d-10 ) exit
       phiG0(:)=phiG(:)

    end do ! i

    deallocate( phiG0 )
    deallocate( phi )

  END SUBROUTINE calc_phiG

  SUBROUTINE check_phiG( nr, dr, phi, phiG )
    implicit none
    integer,intent(IN) :: nr
    real(8),intent(IN) :: dr, phi(nr)
    real(8),intent(OUT) :: phiG(:)
    integer :: i,ig
    real(8) :: r,G,sum0,pi4

    pi4 = 4.0d0*acos(-1.0d0)

    do ig=1,NMGL

       G=sqrt( GG(ig) )

       if ( G <= 1.d-9 ) then

          sum0=0.0d0
          do i=1,nr
             r = i*dr
             sum0 = sum0 + r*r*phi(i)
          end do
          phiG(ig)=pi4*sum0*dr

       else

          sum0=0.0d0
          do i=1,nr
             r = i*dr
             sum0 = sum0 + r*phi(i)*sin(G*r)
          end do
          phiG(ig)=pi4/G*sum0*dr

       end if

    end do ! ig

  END SUBROUTINE check_phiG

  SUBROUTINE make_phiG_by_interpolation( phiGtab, phiGout )
    implicit none
    real(8),intent(IN)  :: phiGtab(:,0:,0:)
    real(8),intent(OUT) :: phiGout(:,0:,0:)
    integer :: ig,i,i0,i1,m,j,k
    real(8) :: G,dG,f,f0,err,err0
    real(8),allocatable :: Gtab(:)
    real(8),parameter :: tol=1.d-10

    dG = Gmax_table/NG_table

    allocate( Gtab(NG_table) ) ; Gtab=0.0d0
    do i=1,NG_table
       Gtab(i) = (i-1)*dG
    end do

    do ig=1,NMGL

       G=sqrt( GG(ig) )

       if ( G < 1.d-9 ) then
          phiGout(ig,:,:) = phiGtab(1,:,:)
       else if ( G > Gmax_table ) then
          phiGout(ig,:,:) = phiGtab(NG_table,:,:)
       else
          i=nint(G/dG)
          do k=1,NumQgrid
          do j=k,NumQgrid
             err0=1.d10
             do m=1,20
                i0=max(i-m,1)
                i1=min(i+m,NG_table)
                call polint(Gtab(i0:i1),phiGtab(i0:i1,j,k),i1-i0+1,G,f,err)
                if ( abs(err) < err0 ) then
                   f0=f
                   err0=abs(err)
                   if ( err0 < tol ) exit
                end if
             end do ! m
             phiG(ig,j,k)=f0
             if ( j /= k ) phiG(ig,k,j)=f0
          end do ! j
          end do ! k
       end if

    end do ! ig

    deallocate( Gtab )

  END SUBROUTINE make_phiG_by_interpolation

  SUBROUTINE calc_phi_b( di, dj, phi )
    implicit none
    real(8),intent(IN)  :: di,dj
    real(8),intent(OUT) :: phi
    integer :: nt,nt0,it,loop
    real(8) :: sum1,phi0,t,dt,dt0,pi
    real(8),parameter :: tol=1.d-5

    pi = acos(-1.0d0)

    nt0 = 8
    dt0 = 0.5d0*pi / nt0

    phi  = 0.0d0
    phi0 = 0.0d0

    nt = nt0/2 ! Interval [0,pi/4] is enough due to symmetry
    dt = dt0

    call Int_RadialPart( 0.25d0*pi, di, dj, sum1 )

    phi = 0.5d0*( 0.0d0 + sum1 )

    do it=1,nt-1
       t = it*dt
       call Int_RadialPart( t, di, dj, sum1 )
       phi = phi + sum1
    end do ! it

    nt = nt/2

    do loop=1,10

       nt = nt*2
       dt = dt*0.5d0

       do it=1,nt
          t = (2*it-1)*dt
          call Int_RadialPart( t, di, dj, sum1 )
          phi = phi + sum1
       end do ! it

       if ( abs(phi*dt-phi0) <= tol ) exit

       phi0 = phi*dt

    end do ! loop

    phi = 2.0d0*(2.0d0/pi**2)*phi*dt

  END SUBROUTINE calc_phi_b

  SUBROUTINE Int_RadialPart( t, di, dj, sum_rad )
    implicit none
    real(8),intent(IN)  :: t,di,dj
    real(8),intent(OUT) :: sum_rad
    integer :: nx0,nx,ix,loop
    real(8) :: sum0,sum1,const
    real(8) :: xmin,xmax,dx,dx0,x,coshx,sinhx,r
    real(8),parameter :: tol=1.d-8

    xmax = 3.0d0
    xmin =-3.0d0
    nx0  = 50
    dx0  = (xmax-xmin)/nx0

    sum0 = 0.0d0
    sum1 = 0.0d0

    const = 0.5d0*acos(-1.0d0)

    nx = nx0
    dx = dx0
    do ix=1,nx-1
       x=xmin+ix*dx
       coshx = 0.5d0*( exp(x) + exp(-x) )
       sinhx = 0.5d0*( exp(x) - exp(-x) )
       r = exp(const*sinhx)
       sum0 = sum0 + phi_b(r,t,di,dj)*coshx*const*r
!       sum0 = sum0 + phi_c(r,t,di,dj)*coshx*const*r
    end do ! ix
    sum1=sum0*dx

    nx = nx/2

    do loop=1,20

       nx = nx*2
       dx = dx*0.5d0
       do ix=1,nx
          x=xmin+(2*ix-1)*dx
          coshx = 0.5d0*( exp(x) + exp(-x) )
          sinhx = 0.5d0*( exp(x) - exp(-x) )
          r = exp(const*sinhx)
          sum0 = sum0 + phi_b(r,t,di,dj)*coshx*const*r
!          sum0 = sum0 + phi_c(r,t,di,dj)*coshx*const*r
       end do ! ix

       if ( abs(sum0*dx-sum1) <= tol ) exit
       sum1=sum0*dx

    end do ! loop

    sum_rad = sum0*dx

  END SUBROUTINE Int_RadialPart

  FUNCTION phi_b( r, theta, di, dj )
    implicit none
    real(8) :: phi_b
    real(8),intent(IN)  :: r, theta, di, dj
    real(8) :: a,b,WWW,TTT !,gamma,pi
    real(8) :: ct,st,sa,sb,ca,cb,hai,haj,hbi,hbj,vai,vaj,vbi,vbj
    real(8) :: a2,b2,cai,cbi,caj,cbj,gdidi_i,gdjdj_i
    real(8),parameter :: pi=3.141592653589793d0
    real(8),parameter :: gamma=1.396263401595464d0
    real(8),parameter :: thr=3.0d0

!    pi = acos(-1.0d0)
!    gamma = 4.0d0*pi/9.0d0

    a = r*cos(theta)
    b = r*sin(theta)

    a2=a*a
    b2=b*b

    gdidi_i = gamma/(di*di)
    gdjdj_i = gamma/(dj*dj)

    cai = a2*gdidi_i
    cbi = b2*gdidi_i
    caj = a2*gdjdj_i
    cbj = b2*gdjdj_i
!    cai = gamma*a2/(di*di)
!    cbi = gamma*b2/(di*di)
!    caj = gamma*a2/(dj*dj)
!    cbj = gamma*b2/(dj*dj)

!    ct = cos(theta)
!    st = sin(theta)

!    a = r*ct
!    b = r*st

!    ca = cos(a)
!    cb = cos(b)
!    sa = 1.0d0 ; if ( a /= 0.0d0 ) sa=sin(a)/a
!    sb = 1.0d0 ; if ( b /= 0.0d0 ) sb=sin(b)/b
!    sa=sin(a)/a
!    sb=sin(b)/b

!    WWW = ( (thr-a2)*sa-thr*ca )*cos(b) &
!        + ( (thr-b2)*ca + (r*r-thr)*sa )*sin(b)/b
!    WWW = ( (thr-a2)*sa-thr*ca )*cb + ( (thr-b2)*ca + (r*r-thr)*sa )*sb
!    WWW = (thr-a2)*cb*sa + (thr-b2)*ca*sb + (r*r-thr)*sa*sb - thr*ca*cb
!    WWW = (3.0d0-a2)*cb*sa + (3.0d0-b2)*ca*sb &
!        + (r*r-3.0d0)*sa*sb - 3.0d0*ca*cb
!    WWW = 2.0d0*( (3.0d0-a2)*cb*sa + (3.0d0-b2)*ca*sb &
!                + (r*r-3.0d0)*sa*sb - 3.0d0*ca*cb )
!    WWW = 2.0d0*( (3.0d0-a*a)*cb*sa + (3.0d0-b*b)*ca*sb &
!                + (r*r-3.0d0)*sa*sb - 3.0d0*ca*cb )

!    if ( di == 0.0d0 ) then
!       hai = 0.0d0
!       hbi = 0.0d0
!    else
!       hai = 1.0d0 - exp( -gamma*(a/di)**2 )
!       hbi = 1.0d0 - exp( -gamma*(b/di)**2 )
!       hai = 1.0d0 - exp( -cai )
!       hbi = 1.0d0 - exp( -cbi )
!    end if
!    if ( dj == 0.0d0 ) then
!       haj = 0.0d0
!       hbj = 0.0d0
!    else
!       haj = 1.0d0 - exp( -gamma*(a/dj)**2 )
!       hbj = 1.0d0 - exp( -gamma*(b/dj)**2 )
!       haj = 1.0d0 - exp( -caj )
!       hbj = 1.0d0 - exp( -cbj )
!    end if

!    if ( a == 0.0d0 ) then
!       vai = di*di/(2.0d0*gamma)
!       vaj = dj*dj/(2.0d0*gamma)
!    else
!       vai = a2/(2.0d0*hai)
!       vaj = a2/(2.0d0*haj)
!       vai = a2/( 2.0d0-2.0d0*exp(-cai) )
!       vaj = a2/( 2.0d0-2.0d0*exp(-caj) )
       vai = 0.5d0*a2/( 1.0d0-exp(-cai) )
       vaj = 0.5d0*a2/( 1.0d0-exp(-caj) )
!       vai = a2/( 1.0d0-exp(-cai) )
!       vaj = a2/( 1.0d0-exp(-caj) )
!    end if
!    if ( b == 0.0d0 ) then
!       vbi = di*di/(2.0d0*gamma)
!       vbj = dj*dj/(2.0d0*gamma)
!    else
!       vbi = b2/(2.0d0*hbi)
!       vbj = b2/(2.0d0*hbj)
!       vbi = b2/( 2.0d0-2.0d0*exp(-cbi) )
!       vbj = b2/( 2.0d0-2.0d0*exp(-cbj) )
       vbi = 0.5d0*b2/( 1.0d0-exp(-cbi) )
       vbj = 0.5d0*b2/( 1.0d0-exp(-cbj) )
!       vbi = b2/( 1.0d0-exp(-cbi) )
!       vbj = b2/( 1.0d0-exp(-cbj) )
!    end if

    TTT = ( 1.0d0/(vai+vbi) + 1.0d0/(vaj+vbj) ) &
               *( 1.0d0/( (vai+vaj)*(vbi+vbj) ) &
                 +1.0d0/( (vai+vbj)*(vaj+vbi) ) )
!    TTT = ( 1.0d0/(vai+vbi) + 1.0d0/(vaj+vbj) ) &
!               *( 1.0d0/( (vai+vaj)*(vbi+vbj) ) &
!                 +1.0d0/( (vai+vbj)*(vaj+vbi) ) )
!    TTT = 0.5d0*( 1.0d0/(vai+vbi) + 1.0d0/(vaj+vbj) ) &
!               *( 1.0d0/( (vai+vaj)*(vbi+vbj) ) &
!                 +1.0d0/( (vai+vbj)*(vaj+vbi) ) )

    ca = cos(a)
    sa = sin(a)/a
    WWW = ( (thr-a2)*sa-thr*ca )*cos(b) &
        + ( (thr-b2)*ca + (r*r-thr)*sa )*sin(b)/b

    phi_b = r*WWW*TTT

  END FUNCTION phi_b

  FUNCTION phi_c( r, theta, di, dj )
    implicit none
    real(8) :: phi_c
    real(8),intent(IN)  :: r, theta, di, dj
    real(8) :: a,b,gamma,pi,WWW,TTT
    real(8) :: ct,st,sa,sb,ca,cb,hai,haj,hbi,hbj,vai,vaj,vbi,vbj

    pi = acos(-1.0d0)
    gamma = 4.0d0*pi/9.0d0

    ct = cos(theta)
    st = sin(theta)

    a = r*ct
    b = r*st

    ca = cos(a)
    cb = cos(b)
    sa = 1.0d0 ; if ( a /= 0.0d0 ) sa=sin(a)/a
    sb = 1.0d0 ; if ( b /= 0.0d0 ) sb=sin(b)/b

    WWW = 2.0d0*( (3.0d0-a*a)*cb*sa + (3.0d0-b*b)*ca*sb &
                + (r*r-3.0d0)*sa*sb - 3.0d0*ca*cb )

    if ( di == 0.0d0 ) then
       hai = 0.0d0
       hbi = 0.0d0
    else
       hai = 1.0d0 - exp( -gamma*(a/di)**2 )
       hbi = 1.0d0 - exp( -gamma*(b/di)**2 )
    end if
    if ( dj == 0.0d0 ) then
       haj = 0.0d0
       hbj = 0.0d0
    else
       haj = 1.0d0 - exp( -gamma*(a/dj)**2 )
       hbj = 1.0d0 - exp( -gamma*(b/dj)**2 )
    end if

    if ( a == 0.0d0 ) then
       vai = di*di/(2.0d0*gamma)
       vaj = dj*dj/(2.0d0*gamma)
    else
       vai = a*a/(2.0d0*hai)
       vaj = a*a/(2.0d0*haj)
    end if
    if ( b == 0.0d0 ) then
       vbi = di*di/(2.0d0*gamma)
       vbj = dj*dj/(2.0d0*gamma)
    else
       vbi = b*b/(2.0d0*hbi)
       vbj = b*b/(2.0d0*hbj)
    end if

    TTT = 0.5d0*( 1.0d0/(vai+vbi) + 1.0d0/(vaj+vbj) ) &
               *( 1.0d0/( (vai+vaj)*(vbi+vbj) ) &
                 +1.0d0/( (vai+vbj)*(vaj+vbi) ) )

    phi_c = r*WWW*TTT

  END FUNCTION phi_c

  SUBROUTINE init1_xc_vdw
    implicit none
    real(8),allocatable :: mat1(:,:),mat2(:,:),work(:)
    integer,allocatable :: ipiv(:)
    real(8) :: h
    integer :: m,n,i,info

    n = NumQgrid
    h = SizeQgrid
    m = n+1

    allocate( bmat(0:n,0:n) ) ; bmat=0.0d0
    allocate( cmat(0:n,0:n) ) ; cmat=0.0d0
    allocate( dmat(0:n,0:n) ) ; dmat=0.0d0

    allocate( mat1(0:n,0:n) ) ; mat1=0.0d0
    allocate( mat2(0:n,0:n) ) ; mat2=0.0d0
    allocate( ipiv(0:n)     ) ; ipiv=0
    allocate( work(64*m)    ) ; work=0.0d0

    mat1(0,0)=1.0d0
    do i=1,n-1
       mat1(i,i-1)=h
       mat1(i,i  )=4.0d0*h
       mat1(i,i+1)=h
    end do
    mat1(n,n)=1.0d0

    do i=1,n-1
       mat2(i,i-1)= 1.0d0
       mat2(i,i  )=-2.0d0
       mat2(i,i+1)= 1.0d0
    end do
    mat2(0:n,0:n)=mat2(0:n,0:n)*3.0d0/h

    call dgetrf( m,m,mat1,m,ipiv,info)
    call dgetri( m,mat1,m,ipiv,work,size(work),info)

    cmat(0:n,0:n) = matmul( mat1(0:n,0:n), mat2(0:n,0:n) )

    mat1=0.0d0
    do i=0,n-1
       mat1(i,i  )=-1.0d0
       mat1(i,i+1)= 1.0d0
    end do

    dmat(0:n,0:n) = matmul( mat1(0:n,0:n), cmat(0:n,0:n) )/(3.0d0*h)

    bmat(0:n,0:n) = mat1(0:n,0:n)/h - cmat(0:n,0:n)*h - dmat(0:n,0:n)*h**2

    deallocate( work )
    deallocate( ipiv )
    deallocate( mat2 )
    deallocate( mat1 )

  END SUBROUTINE init1_xc_vdw

  SUBROUTINE calc_xc_vdw( rgrid, rho, ene, pot )
    implicit none
    type( grid ),intent(IN) :: rgrid
    type( GSArray ),intent(IN) :: rho
    type( xcene ) :: ene
    type( xcpot ) :: pot
    type( gradient16 ) :: grad16
    logical :: disp_sw
    integer :: i,j,m0,m1,mm,nspin
    real(8) :: rho_tmp(2), trho, edx_lda(2), edc_lda,ss_min,ss_max
    real(QP) :: kf,pi,onethr,fouthr,ss,const,c,q0
    real(8) :: q0_min,q0_max,sbuf(2),rbuf(2)
    real(8) :: vdx_lda(2), vdc_lda(2)
    real(8),allocatable :: qsat(:), pol_spline(:,:)
    real(8),allocatable :: ua(:,:), pol_drv(:,:), f(:), g(:)
    complex(8),allocatable :: theta(:,:)

    call write_border( 1, " calc_xc_vdw(start)" )

    pi     = acos(-1.0q0)
    onethr = 1.0q0/3.0q0
    fouthr = 4.0q0/3.0q0
    const  = Zab/9.0q0
    m0     = rgrid%g1%head
    m1     = rgrid%g1%tail
    nspin  = rho%s_range%size_global

    call construct_gradient16( rgrid, rho, grad16 )

    allocate( qsat(m0:m1) ) ; qsat=0.0d0

    Ex_LDA = 0.0d0
    Ec_LDA = 0.0d0

    ss_min = 1.d100
    ss_max =-1.d100

    pot%c%val(:,:)=0.0d0

    do i = pot%xc%g_range%head, pot%xc%g_range%tail

       rho_tmp(1) = rho%val(i,1)
       rho_tmp(2) = rho%val(i,nspin)
       trho       = sum( rho_tmp(1:nspin) )

       if ( trho < zero_density ) cycle

       call calc_edx_lda( nspin, rho_tmp, edx_lda, vdx_lda )
       call calc_edc_lda( nspin, rho_tmp, edc_lda, vdc_lda )

       Ex_LDA = Ex_LDA + sum( rho_tmp(1:nspin)*edx_lda(1:nspin) )
       Ec_LDA = Ec_LDA + trho*edc_lda

       do j = pot%xc%s_range%head, pot%xc%s_range%tail
          pot%c%val(i,j) = edc_lda + trho*vdc_lda(j)
       end do

       kf = ( 3.0q0*pi*pi*trho )**onethr

       if ( trho <= 0.0d0 ) then
          qsat(i)=Qmax
       else
          ss = grad16%gg(i)/(2.0q0*kf*trho)**2
          q0 = -fouthr*pi*(edx_lda(1)+edc_lda) - const*ss*kf
          c=0.0q0
          do j=1,Mmax
             c=c+(q0/Qmax)**j/j
          end do
          qsat(i) = Qmax*( 1.0q0 - exp(-c) )
          ss_min = min( ss, ss_min )
          ss_max = max( ss, ss_max )
       end if

    end do ! i

    sbuf(1) = minval(qsat)
    sbuf(2) = maxval(qsat)
    call MPI_ALLREDUCE(sbuf(1),q0_min,1,MPI_REAL8,MPI_MIN,comm_grid,i)
    call MPI_ALLREDUCE(sbuf(2),q0_max,1,MPI_REAL8,MPI_MAX,comm_grid,i)

    sbuf(1) = ss_min
    sbuf(2) = ss_max
    call MPI_ALLREDUCE(sbuf(1),ss_min,1,MPI_REAL8,MPI_MIN,comm_grid,i)
    call MPI_ALLREDUCE(sbuf(2),ss_max,1,MPI_REAL8,MPI_MAX,comm_grid,i)

    sbuf(1) = Ex_LDA * rgrid%VolumeElement
    sbuf(2) = Ec_LDA * rgrid%VolumeElement
    call MPI_ALLREDUCE(sbuf,rbuf,2,MPI_REAL8,MPI_SUM,comm_grid,i)
    Ex_LDA = rbuf(1)
    Ec_LDA = rbuf(2)

!    call check_disp_switch( disp_sw, 0 )
!    if ( disp_sw ) then
!       write(*,*) "ss(min,max)=",ss_min,ss_max
!       write(*,*) "q0(min,max)=",q0_min,q0_max
!       write(*,*) "Ex_LDA =",Ex_LDA
!       write(*,*) "Ec_LDA =",Ec_LDA
!       write(*,*) "Exc_LDA=",Ex_LDA+Ec_LDA
!    end if

! ---

    allocate( pol_spline(m0:m1,0:NumQgrid) ) ; pol_spline=0.0d0

    call calc_pol_spline( m0, m1, qsat, pol_spline )

! ---

    mm = max( m1-m0+1, MG_1-MG_0+1 )

    allocate( theta(mm,0:NumQgrid) ) ; theta=(0.0d0,0.0d0)

    do j=0,NumQgrid
       do i=m0,m1
          theta(i-m0+1,j) = sum( rho%val(i,:) )*pol_spline(i,j)
       end do
    end do

    call fft_theta( rgrid, theta )

! ---

    call calc_vdw_energy( rgrid, theta, Ec_vdw )

    ene%Ec = Ec_LDA + Ec_VDW

! ---

    allocate( ua(m0:m1,0:NumQgrid) ) ; ua=0.0d0

    call calc_ua( rgrid, theta, ua )

! ---

    deallocate( theta )

! ---

    allocate( pol_drv(m0:m1,0:NumQgrid) ) ; pol_drv=0.0d0

    call calc_pol_drv( m0, m1, qsat, pol_spline, pol_drv )

! ---

    allocate( f(m0:m1) ) ; f=0.0d0
    allocate( g(m0:m1) ) ; g=0.0d0

    do j=0,NumQgrid
       f(:) = f(:) + ua(:,j)*pol_drv(:,j)
    end do

    call calc_dqdn( rgrid, rho, grad16, f, g )

! ---

    do j=0,NumQgrid
       g(:) = g(:) + ua(:,j)*pol_spline(:,j)
    end do

    c = 1.0d0/rgrid%g1%size_global

    pot%c%val(:,1) = pot%c%val(:,1) + c*g(:)

! ---

    deallocate( g )
    deallocate( f )
    deallocate( pol_drv )
    deallocate( ua )
    deallocate( pol_spline )
    deallocate( qsat )

    call destruct_gradient16( grad16 )

    call write_border( 1, " calc_xc_vdw(end)" )

  END SUBROUTINE calc_xc_vdw

  SUBROUTINE calc_edx_lda( nn, rho, edx_lda, vdx_lda )
    implicit none
    integer,intent(IN)  :: nn
    real(8),intent(IN)  :: rho(nn)
    real(8),optional,intent(OUT) :: edx_lda(nn)
    real(8),optional,intent(OUT) :: vdx_lda(nn)
    real(8) :: onethr, twothr, thrfou, thrpi, pi, factor
    integer :: s

    pi     = acos(-1.0d0)
    onethr = 1.0d0/3.0d0
    twothr = 2.0d0/3.0d0
    thrfou = 3.0d0/4.0d0
    thrpi  = 3.0d0/pi

    factor = 1.0d0
    if ( nn == 2 ) factor = 2.0d0**onethr

    if ( present(edx_lda) ) then
       do s=1,nn
          edx_lda(s) = -factor*thrfou*( thrPi*rho(s) )**onethr
       end do
    end if

    if ( present(vdx_lda) ) then
       do s=1,nn
          if ( rho(s) > 0.0d0 ) then
             vdx_lda(s) = -pi*thrfou/(3.0d0*pi*pi*rho(s))**twothr
          else
             vdx_lda(s)=0.0d0
          end if
       end do
    end if

  END SUBROUTINE calc_edx_lda

  SUBROUTINE calc_edc_lda( nn, rho, edc_lda, vdc_lda )
    implicit none
    integer,intent(IN)  :: nn
    real(8),intent(IN)  :: rho(nn)
    real(8),optional,intent(OUT) :: edc_lda
    real(8),optional,intent(OUT) :: vdc_lda(nn)
    real(8),parameter :: A00  =0.031091d0,A01  =0.015545d0,A02  =0.016887d0
    real(8),parameter :: alp10=0.21370d0 ,alp11=0.20548d0 ,alp12=0.11125d0
    real(8),parameter :: bt10 =7.5957d0  ,bt11 =1.41189d1 ,bt12 =1.0357d1
    real(8),parameter :: bt20 =3.5876d0  ,bt21 =6.1977d0  ,bt22 =3.6231d0
    real(8),parameter :: bt30 =1.6382d0  ,bt31 =3.3662d0  ,bt32 =0.88026d0
    real(8),parameter :: bt40 =0.49294d0 ,bt41 =0.62517d0 ,bt42 =0.49671d0
    real(8) :: one,two,fouthr,pi,const1,const2,onethr,ThrFouPi
    real(8) :: trho,rhoa,rhob,zeta,fz,kf,rs,rssq,ec_U,ec_P,alpc
    real(8) :: deU_dn, deP_dn, dac_dn, dfz_dz, dz_dn(2), dec_dz
    integer :: s

    trho = sum( rho(1:nn) )
    rhoa = rho(1)
    rhob = rho(nn)

    if ( present(edc_lda) ) edc_lda=0.0d0
    if ( present(vdc_lda) ) vdc_lda=0.0d0

    if ( trho <= zero_density ) return

    one      = 1.0d0
    two      = 2.0d0
    fouthr   = 4.0d0/3.0d0
    onethr   = 1.0d0/3.0d0
    pi       = acos(-1.0d0)
    ThrFouPi = 3.0d0/(4.0d0*pi)
    const1   = two**fouthr-two
    const2   = 9.0d0/4.0d0*(two**onethr-one)

    zeta = ( rhoa - rhob )/trho
    fz   = ( (one+zeta)**fouthr + (one-zeta)**fouthr - two )*const1
    kf   = ( 3.0d0*Pi*Pi*trho )**onethr
    rs   = ( ThrFouPi/trho )**onethr
    rssq = sqrt(rs)
    ec_U = -two*A00*( one + alp10*rs ) &
         *log( one + one/( two*A00 &
         *(bt10*rssq + bt20*rs + bt30*rs*rssq + bt40*rs*rs) ) )
    ec_P = -two*A01*( one + alp11*rs ) &
         *log( one + one/( two*A01 &
         *(bt11*rssq + bt21*rs + bt31*rs*rssq + bt41*rs*rs) ) )
    alpc = -two*A02*( one + alp12*rs ) &
         *log( one + one/( two*A02 &
         *(bt12*rssq + bt22*rs + bt32*rs*rssq + bt42*rs*rs) ) )

    if ( present(edc_lda) ) then
       edc_lda = ec_U - alpc*fz*const2*(one-zeta**4) &
            + ( ec_P - ec_U )*fz*zeta**4
    end if

    if ( present(vdc_lda) ) then

       deU_dn = -4.0d0*Pi/9.0d0*rs**4*alp10*ec_U/( 1.0d0 + alp10*rs ) &
            -4.0d0*Pi/9.0d0*rs*rs*( 1.0d0 + alp10*rs ) &
            *( 0.5d0*bt10*sqrt(rs) + bt20*rs + 1.5d0*bt30*rs*sqrt(rs) &
            +2.0d0*bt40*rs*rs ) &
            /( bt10 + bt20*sqrt(rs) + bt30*rs &
            + bt40*rs*sqrt(rs))**2 &
            * exp( ec_U/(2.0d0*A00*(1.0d0+alp10*rs)) )

       deP_dn = -4.0d0*Pi/9.0d0*rs**4*alp11*ec_P/( 1.0d0 + alp11*rs ) &
            -4.0d0*Pi/9.0d0*rs*rs*( 1.0d0 + alp11*rs ) &
            *( 0.5d0*bt11*sqrt(rs) + bt21*rs + 1.5d0*bt31*rs*sqrt(rs) &
            +2.0d0*bt41*rs*rs ) &
            /( bt11 + bt21*sqrt(rs) + bt31*rs &
            + bt41*rs*sqrt(rs))**2 &
            * exp( ec_P/(2.0d0*A01*(1.0d0+alp11*rs)) )

       dac_dn = -4.0d0*Pi/9.0d0*rs**4*alp12*alpc/( 1.0d0 + alp12*rs ) &
            -4.0d0*Pi/9.0d0*rs*rs*( 1.0d0 + alp12*rs ) &
            *( 0.5d0*bt12*sqrt(rs) + bt22*rs + 1.5d0*bt32*rs*sqrt(rs) &
            +2.0d0*bt42*rs*rs ) &
            /( bt12 + bt22*sqrt(rs) + bt32*rs &
            + bt42*rs*sqrt(rs))**2 &
            * exp( alpc/(2.0d0*A02*(1.d0+alp12*rs)) )

       dfz_dz = fouthr*( (one+zeta)**onethr - (one-zeta)**onethr )*const1

       dz_dn(1)  = 2.0d0*rhob/trho
       dz_dn(nn) =-2.0d0*rhoa/trho

       dec_dz = -alpc*dfz_dz*const2*(one-zeta**4) &
               + 4.0d0*alpc*fz*const2*zeta**3 &
               +(ec_P-ec_U)*dfz_dz*zeta**4 &
               +(ec_P-ec_U)*fz*4.0d0*zeta**3

       do s=1,nn
          vdc_lda(s) = deU_dn - dac_dn*fz*const2*(one-zeta**4) &
               +(deP_dn-deU_dn)*fz*zeta**4 + dec_dz*dz_dn(s)
       end do

    end if

  END SUBROUTINE calc_edc_lda

  SUBROUTINE calc_pol_spline( m0, m1, q0, pol )
    implicit none
    integer,intent(IN)  :: m0,m1
    real(8),intent(IN)  :: q0(m0:m1)
    real(8),intent(OUT) :: pol(m0:m1,0:NumQgrid)
    integer :: i,j,k
    real(8) :: q

    pol(:,:)=0.0d0

    do k=m0,m1

       do i=0,NumQgrid
          if ( q0(k) < Qgrid(i) ) exit
       end do
       i=i-1

       q=q0(k)-Qgrid(i)
       do j=0,NumQgrid
          pol(k,j) = bmat(i,j)*q + cmat(i,j)*q**2 + dmat(i,j)*q**3
       end do
       pol(k,i) = pol(k,i) + 1.0d0

! ---------------------------- Wu & Gygi ---
       if ( q0(k) /= 0.0d0 ) then
          do j=0,NumQgrid
             pol(k,j) = Qgrid(j)*pol(k,j)/q0(k)
          end do
       else
          pol(k,:) = 0.0d0
       end if
! ------------------------------------------

    end do ! k

  END SUBROUTINE calc_pol_spline

  SUBROUTINE fft_theta( rgrid, theta )
    implicit none
    type( grid ),intent(IN) :: rgrid
    complex(8),intent(INOUT) :: theta(:,0:)
    integer :: ML,ML1,ML2,ML3,mm,i1,i2,i3,i,j,info
    complex(8),allocatable :: work0(:,:,:), work1(:,:,:)
    complex(8),parameter :: zero=(0.0d0,0.0d0)

    ML  = rgrid%g1%size_global
    ML1 = rgrid%g3%x%size_global
    ML2 = rgrid%g3%y%size_global
    ML3 = rgrid%g3%z%size_global

    allocate( work0(0:ML1-1,0:ML2-1,0:ML3-1) ) ; work0=zero
    allocate( work1(0:ML1-1,0:ML2-1,0:ML3-1) ) ; work1=zero

    call init_fft

    call construct_ggrid(1)

    do j=0,NumQgrid

       work1(:,:,:)=zero
       i=0
       do i3=rgrid%g3%z%head,rgrid%g3%z%tail
       do i2=rgrid%g3%y%head,rgrid%g3%y%tail
       do i1=rgrid%g3%x%head,rgrid%g3%x%tail
          i=i+1
          work1(i1,i2,i3) = theta(i,j)
       end do ! i1
       end do ! i2
       end do ! i3

       call MPI_ALLREDUCE( work1, work0, size(work0) &
            , MPI_COMPLEX16, MPI_SUM, comm_grid, info )

       call forward_fft( work0, work1 )

       theta(:,j) = zero
       do i=MG_0,MG_1
          i1=mod( LLG(1,i)+ML1, ML1 )
          i2=mod( LLG(2,i)+ML2, ML2 )
          i3=mod( LLG(3,i)+ML3, ML3 )
          theta(i-MG_0+1,j) = work0(i1,i2,i3)
       end do

    end do ! j

    call destruct_ggrid

    call finalize_fft

    deallocate( work1 )
    deallocate( work0 )

  END SUBROUTINE fft_theta

  SUBROUTINE calc_vdw_energy( rgrid, theta, Ec )
    implicit none
    complex(8),intent(IN) :: theta(:,0:)
    type( grid ),intent(IN) :: rgrid
    real(8),intent(INOUT) :: Ec
    integer :: a,b,i,j,info
    complex(8) :: z
    real(8) :: c,E0
    logical :: disp_sw

    z=(0.0d0,0.0d0)

    do b=0,NumQgrid
    do a=b,NumQgrid

       c=1.0d0 ; if ( a /= b ) c=2.0d0

       do i=MG_0,MG_1
          j=i-MG_0+1
          z = z + c*theta(j,a)*conjg(theta(j,b))*phiG(MGL(i),a,b)
       end do ! i

    end do ! a
    end do ! b

    E0 = 0.5d0*z * rgrid%g1%size_global * rgrid%VolumeElement
    call MPI_ALLREDUCE(E0,Ec,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,info)

!    call check_disp_switch( disp_sw, 0 )
!    if ( disp_sw ) then
!       write(*,*) "Ec_vdw=",Ec
!       write(*,*) "Ec_vdw+Ec_LDA=",Ec+Ec_LDA
!    end if

  END SUBROUTINE calc_vdw_energy

  SUBROUTINE calc_ua( rgrid, theta, ua )
    implicit none
    complex(8),intent(IN) :: theta(:,0:)
    type( grid ),intent(IN) :: rgrid
    real(8),intent(OUT)   :: ua(:,0:)
    integer :: a,b,i,j,i1,i2,i3
    complex(8),parameter :: zero=(0.0d0,0.0d0)
    complex(8),allocatable :: work0(:,:,:),work1(:,:,:)
    integer :: ML1,ML2,ML3,ML,info

    ML  = rgrid%g1%size_global
    ML1 = rgrid%g3%x%size_global
    ML2 = rgrid%g3%y%size_global
    ML3 = rgrid%g3%z%size_global

    allocate( work0(0:ML1-1,0:ML2-1,0:ML3-1) ) ; work0=zero
    allocate( work1(0:ML1-1,0:ML2-1,0:ML3-1) ) ; work1=zero

    call init_fft

    call construct_ggrid(1)

    ua(:,:) = 0.0d0

    do a=1,NumQgrid

       work1(:,:,:)=zero

       do b=1,NumQgrid

          do i=MG_0,MG_1
             j=i-MG_0+1
             i1 = mod( ML1+LLG(1,i), ML1 )
             i2 = mod( ML2+LLG(2,i), ML2 )
             i3 = mod( ML3+LLG(3,i), ML3 )
             work1(i1,i2,i3) = work1(i1,i2,i3) + phiG(MGL(i),a,b)*theta(j,b)
          end do ! i

       end do ! b

       call MPI_ALLREDUCE(work1,work0,size(work0),MPI_COMPLEX16 &
            ,MPI_SUM,MPI_COMM_WORLD,info)

       call backward_fft( work0, work1 )

       i=0
       do i3=rgrid%g3%z%head,rgrid%g3%z%tail
       do i2=rgrid%g3%y%head,rgrid%g3%y%tail
       do i1=rgrid%g3%x%head,rgrid%g3%x%tail
          i=i+1
          ua(i,a) = work0(i1,i2,i3)
       end do
       end do
       end do

    end do ! a

    call destruct_ggrid

    call finalize_fft

    deallocate( work1 )
    deallocate( work0 )

  END SUBROUTINE calc_ua

  SUBROUTINE calc_pol_drv( m0, m1, q0, pol, pol_drv )
    implicit none
    integer,intent(IN)  :: m0,m1
    real(8),intent(IN)  :: q0(m0:m1),pol(m0:m1,0:NumQgrid)
    real(8),intent(OUT) :: pol_drv(m0:m1,0:NumQgrid)
    integer :: i,j,k
    real(8) :: q

    pol_drv(:,:)=0.0d0

    do k=m0,m1

       do i=0,NumQgrid
          if ( q0(k) < Qgrid(i) ) exit
       end do
       i=i-1

       q=q0(k)-Qgrid(i)
       do j=0,NumQgrid
          pol_drv(k,j) = bmat(i,j) + 2.0d0*cmat(i,j)*q + 3.0d0*dmat(i,j)*q*q
       end do

! ---------------------------- Wu & Gygi ---
       if ( q0(k) /= 0.0 ) then
          do j=0,NumQgrid
             pol_drv(k,j)=Qgrid(j)*(-pol(k,j)+q0(k)*pol_drv(k,j))/q0(k)**2
          end do
       else
          pol_drv(k,:)=0.0d0
       end if
! ------------------------------------------

    end do ! k

  END SUBROUTINE calc_pol_drv

  SUBROUTINE calc_dqdn( rgrid, rho, grad16, fin, fou )
    implicit none
    type( grid ) :: rgrid
    type( GSArray ) :: rho
    type( gradient16 ) :: grad16
    real(8),intent(IN)  :: fin(:)
    real(8),intent(OUT) :: fou(:)
    type( fd ) :: nabla
    type( lattice ) :: aa, bb
    real(8) :: trho,rho_tmp(2),vdx_lda(2),vdc_lda(2)
    real(8) :: edx_lda(2),edc_lda
    real(QP) :: pi,FouPiThr,twothr,onethr,sevthr,kf,ss,cm,q0,c,d,const
    real(QP) :: qtrho,b(3,3)
    real(QP),allocatable :: rtmp0(:),rtmp(:,:),ftmp(:)
    integer :: m0,m1,m,ML1,ML2,ML3,Md,i,i1,i2,i3,j,j1,j2,j3,k1,k2,k3,info
    integer :: nspin
    integer,allocatable :: LLL(:,:,:)

    pi       = acos(-1.0q0)
    onethr   = 1.0q0/3.0q0
    twothr   = 2.0q0/3.0q0
    FouPiThr = 4.0q0*pi/3.0q0
    sevthr   = 7.0q0/3.0q0
    const    = Zab/9.0q0

    m0 = rgrid%g1%head
    m1 = rgrid%g1%tail

    ML1 = rgrid%g3%x%size_global
    ML2 = rgrid%g3%y%size_global
    ML3 = rgrid%g3%z%size_global

    nspin = rho%s_range%size_global

    allocate( rtmp(rgrid%g1%size_global,3) ) ; rtmp =0.0q0
    allocate( rtmp0(m0:m1)                 ) ; rtmp0=0.0q0
    allocate( ftmp(m0:m1)                  ) ; ftmp =0.0q0

! ---

    do i=m0,m1

       rho_tmp(1) = rho%val(i,1)
       rho_tmp(2) = rho%val(i,nspin)
       trho       = sum( rho_tmp(1:nspin) )
       qtrho      = trho

       if ( trho <= 0.0d0 ) cycle
!       if ( trho <= zero_density ) cycle

       call calc_edx_lda( nspin, rho_tmp, edx_lda, vdx_lda )
       call calc_edc_lda( nspin, rho_tmp, edc_lda, vdc_lda )

       kf = (3.0q0*pi*pi*qtrho)**onethr
       ss = grad16%gg(i)/(2.0q0*kf*qtrho)**2
       q0 = -FouPithr*(edx_lda(1)+edc_lda) - const*ss*kf

       c=1.0q0
       do j=1,Mmax-1
          c=c+(q0/Qmax)**j
       end do
       d=0.0q0
       do j=1,Mmax
          d=d+(q0/Qmax)**j/j
       end do
       c=c*exp(-d)

       ftmp(i) = c*( -FouPiThr*( vdx_lda(1) + vdc_lda(1) )*trho &
                     +const*sevthr*kf*ss )*fin(i-m0+1)

       rtmp0(i) = c*( -const/(2.0q0*qtrho*kf) )*fin(i-m0+1)

    end do ! i

! ---

    rtmp(m0:m1,1) = rtmp0(m0:m1) * grad16%gx(m0:m1)
    rtmp(m0:m1,2) = rtmp0(m0:m1) * grad16%gy(m0:m1)
    rtmp(m0:m1,3) = rtmp0(m0:m1) * grad16%gz(m0:m1)

    deallocate( rtmp0 )

    do i=1,3
       call MPI_ALLGATHERV( rtmp(m0,i), m1-m0+1, MPI_REAL16, rtmp(1,i) &
            , ir_grid, id_grid, MPI_REAL16, comm_grid, info )
    end do

    allocate( LLL(0:ML1-1,0:ML2-1,0:ML3-1) ) ; LLL=0

    i=m0-1
    do i3=rgrid%g3%z%head,rgrid%g3%z%tail
    do i2=rgrid%g3%y%head,rgrid%g3%y%tail
    do i1=rgrid%g3%x%head,rgrid%g3%x%tail
       i=i+1
       LLL(i1,i2,i3) = i
    end do
    end do
    end do
    call MPI_ALLREDUCE( MPI_IN_PLACE, LLL, size(LLL), MPI_INTEGER &
         ,MPI_SUM, comm_grid, info )

    call construct_nabla_fd( nabla )
    Md = nabla%md

    call get_aa_lattice( aa )
    call get_reciprocal_lattice( aa, bb )
    b(1:3,1)=aa%Length(1)*bb%LatticeVector(1:3,1)/( 2*pi*rgrid%spacing(1) )
    b(1:3,2)=aa%Length(2)*bb%LatticeVector(1:3,2)/( 2*pi*rgrid%spacing(2) )
    b(1:3,3)=aa%Length(3)*bb%LatticeVector(1:3,3)/( 2*pi*rgrid%spacing(3) )

    do i3=0,ML3-1
    do i2=0,ML2-1
    do i1=0,ML1-1
       i=LLL(i1,i2,i3)
       do m=-Md,Md
          cm=nabla%coef(m)*sign(1,m)
          j1=i1+m
          k1=j1/ML1 ; if ( j1 < 0 ) k1=(j1+1)/ML1-1
          j1=j1-k1*ML1
          j =LLL(j1,i2,i3)
          if ( m0 <= j .and. j <= m1 ) then
             ftmp(j) = ftmp(j) + cm*( rtmp(i,1)*b(1,1) &
                                     +rtmp(i,2)*b(2,1) &
                                     +rtmp(i,3)*b(3,1) )
          end if
          j2=i2+m
          k2=j2/ML2 ; if ( j2 < 0 ) k2=(j2+1)/ML2-1
          j2=j2-k2*ML2
          j =LLL(i1,j2,i3)
          if ( m0 <= j .and. j <= m1 ) then
             ftmp(j) = ftmp(j) + cm*( rtmp(i,1)*b(1,2) &
                                     +rtmp(i,2)*b(2,2) &
                                     +rtmp(i,3)*b(3,2) )
          end if
          j3=i3+m
          k3=j3/ML3 ; if ( j3 < 0 ) k3=(j3+1)/ML3-1
          j3=j3-k3*ML3
          j =LLL(i1,i2,j3)
          if ( m0 <= j .and. j <= m1 ) then
             ftmp(j) = ftmp(j) + cm*( rtmp(i,1)*b(1,3) &
                                     +rtmp(i,2)*b(2,3) &
                                     +rtmp(i,3)*b(3,3) )
          end if
       end do ! m
    end do ! i1
    end do ! i2
    end do ! i3

    fou(1:m1-m0+1) = ftmp(m0:m1)

! ---

    deallocate( ftmp )
    deallocate( LLL  )
    deallocate( rtmp )

    call destruct_nabla_fd( nabla )

  END SUBROUTINE calc_dqdn

END MODULE xc_vdw_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE BD_CPENTA ( CARD, LARGE_FLD_INP, NUM_GRD )

! Processes CPENTA Bulk Data Cards
!  1) Sets ETYPE for this element type
!  2) Calls subr ELEPRO to read element ID, property ID and connection data into array EDAT

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  FATAL_ERR, JCARD_LEN, NCPENTA6, NCPENTA15, NEDAT, NELE, BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_CPENTA_BEGEND
      USE MODEL_STUF, ONLY            :  ETYPE

      USE BD_CPENTA_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BD_CPENTA'
      CHARACTER(LEN=*), INTENT(INOUT) :: CARD              ! A Bulk Data card
      CHARACTER(LEN=*), INTENT(IN)    :: LARGE_FLD_INP     ! If 'Y', CARD is large field format
      CHARACTER(LEN(CARD))            :: CHILD             ! "Child" card read in subr NEXTC, called herein
      CHARACTER(LEN=JCARD_LEN)        :: JCARD(10)         ! The 10 fields of characters making up CARD
      CHARACTER(LEN(JCARD))           :: ID                ! Character value of element ID (field 2 of parent card)
      CHARACTER(LEN(JCARD))           :: JCARD_EDAT(10)    ! JCARD values sent to subr ELEPRO
      CHARACTER(LEN(JCARD))           :: NAME              ! JCARD(1) from parent entry

      INTEGER(LONG), INTENT(OUT)      :: NUM_GRD           ! Number of GRID's + SPOINT's for the elem
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: ICONT     = 0     ! Indicator of whether a cont card exists. Output from subr NEXTC
      INTEGER(LONG)                   :: IERR      = 0     ! Error indicator returned from subr NEXTC called herein
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_CPENTA_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! CPENTA element Bulk Data Card routine

!   FIELD   ITEM                   ARRAY ELEMENT
!   -----   ------------   ---------------------------------

! Mandatory parent card

!    1      Element type   ETYPE(nele) = 'PENTA6' or 'PENTA15'
!    2      Element ID     EDAT(nedat+1)
!    3      Property ID    EDAT(nedat+1)
!    4-9    Grids 1-6      EDAT(nedat+3) thru EDAT(nedat+8)

! Possible 2nd card if this is a 15 node PENTA:

!    2-9    Grids 7-14     EDAT(nedat+9) thru EDAT(nedat+16) if this is a PENTA 15 node element

! Possible 3rd card if this is a 15 node PENTA:

!    2      Grid  15       EDAT(nedat+17) if this is a PENTA 15 node element

! Make JCARD from CARD

      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
      NAME = JCARD(1)
      ID   = JCARD(2)

! Set JCARD_EDAT to JCARD

      DO I=1,10
         JCARD_EDAT(I) = JCARD(I)
      ENDDO

! Read and check data

      CALL ELEPRO ( 'Y', JCARD_EDAT, 8, 8, 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y' )

      CALL BD_IMBEDDED_BLANK   ( JCARD,2,3,4,5,6,7,8,9 )   ! Make sure that there are no imbedded blanks in fields 2-9
      CALL CRDERR ( CARD )                                 ! CRDERR prints errors found when reading fields

! Check for continuation (if so must be a PENTA15 which then requires 3rd card also)

      ETYPE(NELE)(1:) = ' '
      IF (LARGE_FLD_INP == 'N') THEN
         CALL NEXTC  ( CARD, ICONT, IERR )
      ELSE
         CALL NEXTC2 ( CARD, ICONT, IERR, CHILD )
         CARD = CHILD
      ENDIF
      CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
      IF (ICONT == 0) THEN                                 ! Assume PENTA6 since there was no 2nd card

         NCPENTA6    = NCPENTA6 + 1                        ! Default values unless there is a non-blank continuation
         ETYPE(NELE) = 'PENTA6   '
         NUM_GRD     = 6

      ELSE

         IF (CARD(1:) /= ' ') THEN                         ! Only process continuation entries if 1st one is not totally blank
                                                           ! Assume PENTA15 since there was a 2nd card
            NCPENTA15   = NCPENTA15 + 1
            ETYPE(NELE) = 'PENTA15  '
            NUM_GRD     = 15

            DO I=1,10
               JCARD_EDAT(I) = JCARD(I)
            ENDDO

            CALL ELEPRO ( 'N', JCARD_EDAT, 8, 8, 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y', 'Y' )
            CALL BD_IMBEDDED_BLANK ( JCARD,2,3,4,5,6,7,8,9 )
            CALL CRDERR ( CARD )

            IF (LARGE_FLD_INP == 'N') THEN
               CALL NEXTC  ( CARD, ICONT, IERR )           ! Required 3rd card for PENTA15
            ELSE
               CALL NEXTC2 ( CARD, ICONT, IERR, CHILD )
               CARD = CHILD
            ENDIF
            CALL MKJCARD ( SUBR_NAME, CARD, JCARD )
            IF (ICONT == 1) THEN

               DO I=1,10
                  JCARD_EDAT(I) = JCARD(I)
               ENDDO

               CALL ELEPRO ( 'N', JCARD_EDAT, 1, 1, 'Y', 'N', 'N', 'N', 'N', 'N', 'N', 'N' )

               CALL BD_IMBEDDED_BLANK( JCARD,2,0,0,0,0,0,0,0 )
               CALL CARD_FLDS_NOT_BLANK(JCARD,0,3,4,5,6,7,8,9)
               CALL CRDERR ( CARD )

            ELSE

               FATAL_ERR = FATAL_ERR + 1
               WRITE(ERR,1136) NAME, ID
               WRITE(F06,1136) NAME, ID

            ENDIF

         ELSE                                              ! Since no cont data was read, default to PENTA6

            NCPENTA6    = NCPENTA6 + 1
            ETYPE(NELE) = 'PENTA6   '
            NUM_GRD     = 6
            ENDIF

      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1136 FORMAT(' *ERROR  1136: REQUIRED CONTINUATION FOR ',A,' ID = ',A,' MISSING')

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

      END SUBROUTINE BD_CPENTA


text:
!========================================
! Marc 20/5/03
! Given the netcdf ID, ncid, this file will read variable from it
! File must first be opened using subroutine open_file
!========================================
      subroutine get1d_data_nc(ncid, varname, dim1
     &     , arrayout,ifail)
      implicit none

!----------------------------------------
! For precision
!----------------------------------------
      integer :: realkind
!----------------------------------------
! Define variables in the subroutine header
!----------------------------------------
      integer, intent(in)      :: ncid      ! netCDF dataset ID
      character(*), intent(in) :: varname   ! name of variable to collect
      integer, intent(in)      :: dim1      ! Size of array in 1st dimension
      real, intent(out) :: arrayout(dim1) ! the output array
!----------------------------------------
! Define other netcdf variables
!----------------------------------------
      integer :: status                     ! return code
     &     , varid                          ! ID of variable
!     &     , dim1id                         ! ID for 1st dimension
!----------------------------------------
! local variables
!----------------------------------------
      integer :: ndims                      ! dims. for netcdf variable
     &     , dimid                          ! dimension IDs
     &     , dim1nc                         ! 1st dim. of netcdf var.
!----------------------------------------
! Include files
!----------------------------------------
      include 'netcdf.inc'
      integer ifail

      realkind=kind(arrayout)

!----------------------------------------
! Find the variable ID
!----------------------------------------
      ifail=0
      status=nf_inq_varid(ncid, varname, varid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find ', varname
         ifail=1
         return
      endif
!----------------------------------------
! get information on variable
!----------------------------------------
      status=nf_inq_varndims(ncid, varid, ndims)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find dims. of variable'
         stop
      endif
      if (ndims.ne.1) then
         write(6, *) 'ERROR: variable has ', ndims,
     &        ' dimensions and we expect 1'
      endif
      status=nf_inq_vardimid(ncid, varid, dimid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find dimension ID'
         stop
      endif
!----------------------------------------
! check that dimensions match
!----------------------------------------
      status=nf_inq_dimlen(ncid, dimid, dim1nc)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: Could not get 1st dimension from '
     &        //'netcdf file'
         stop
      endif
      if (dim1nc.ne.dim1) then
         write(6, *) 'ERROR: 1st dimension of variable in model '
     &        //'and netcdf file do not match'
         write(6, *) 'model and netcdf dims are ', dim1, ' and '
     &        , dim1nc
         stop
      endif
!----------------------------------------
! get variable
!----------------------------------------
      if (realkind.eq.4) then
        status=nf_get_var_real(ncid, varid, arrayout)
      else if (realkind.eq.8) then
        status=nf_get_var_double(ncid, varid, arrayout)
      else
        print*,'precision problem in get1d_data_nc'
        stop
      endif
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: getting variable'
         stop
      endif

      end
!
!     Integer version of above code. Any changes in above should
!     also be reflected here
!
      subroutine get1di_data_nc(ncid, varname, dim1
     &     , arrayout,ifail)
      implicit none
!----------------------------------------
! Define variables in the subroutine header
!----------------------------------------
      integer, intent(in)      :: ncid      ! netCDF dataset ID
      character(*), intent(in) :: varname   ! name of variable to collect
      integer, intent(in)      :: dim1      ! Size of array in 1st dimension
      integer, intent(out) :: arrayout(dim1) ! the output array
!----------------------------------------
! Define other netcdf variables
!----------------------------------------
      integer :: status                     ! return code
     &     , varid                          ! ID of variable
!     &     , dim1id                         ! ID for 1st dimension
!----------------------------------------
! local variables
!----------------------------------------
      integer :: ndims                      ! dims. for netcdf variable
     &     , dimid                          ! dimension IDs
     &     , dim1nc                         ! 1st dim. of netcdf var.
!----------------------------------------
! Include files
!----------------------------------------
      include 'netcdf.inc'
      integer ifail
!----------------------------------------
! Find the variable ID
!----------------------------------------
      ifail=0
      status=nf_inq_varid(ncid, varname, varid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find ', varname
         ifail=1
         return
      endif
!----------------------------------------
! get information on variable
!----------------------------------------
      status=nf_inq_varndims(ncid, varid, ndims)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find dims. of variable'
         stop
      endif
      if (ndims.ne.1) then
         write(6, *) 'ERROR: variable has ', ndims,
     &        ' dimensions and we expect 1'
      endif
      status=nf_inq_vardimid(ncid, varid, dimid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find dimension ID'
         stop
      endif
!----------------------------------------
! check that dimensions match
!----------------------------------------
      status=nf_inq_dimlen(ncid, dimid, dim1nc)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: Could not get 1st dimension from '
     &        //'netcdf file'
         stop
      endif
      if (dim1nc.ne.dim1) then
         write(6, *) 'ERROR: 1st dimension of variable in model '
     &        //'and netcdf file do not match'
         write(6, *) 'model and netcdf dims are ', dim1, ' and '
     &        , dim1nc
         stop
      endif
!----------------------------------------
! get variable
!----------------------------------------
      status=nf_get_var_int(ncid, varid, arrayout)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: getting variable'
         stop
      endif

      end

ext: f


text:
El Macero Country Club was founded in 1961 by C. Bruce Mace as part of the El Macero property development. In 1966, the members purchased the club from the Mace family and it became a proprietary position facility as a California NonProfit Corporation. El Macero Country Club is located near Davis, about 10 miles west of Sacramento, in Yolo County. The club, and its golf course, meander through the El Macero development and encompass about 168 acres within beautiful residential surroundings. The entrance to the club is off Mace Boulevard, approximately a mile south of Interstate 80.

This is the Secret clubhouse for everyone who golfs at El Macero.  In addition to the 18hole championship golf course, the club has a heated swimming pool, two recently resurfaced tennis courts, a formal and casual dining area and social events through out the year.  The Club offers Proprietary Golf, Junior Golf and Social NonGolf Memberships. All memberships require an initiation fee and monthly dues and food and beverage minimums vary based upon which membership you are interested in.

Many often assume that you have to be a resident of El Maceroto join the El Macero Country Club; however that is not the case.  Anyone who has enough $$ can join El Macero Country Club, and members come from all over the Sacramento region.

El Macero Country Club is also open to outside events including weddings, meetings and tournaments.  There is an event coordinator on site who is able to assist with all of your planning needs.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
MODULE MODULE_COM

CONTAINS

  FUNCTION GEOMETRICAL_CENTER(coors, atoms_indices_serialized, atoms_series_starts, frame_indices,&
                              n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frame_indices) RESULT(center)

    IMPLICIT NONE

    INTEGER,INTENT(IN)::n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frame_indices
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms,3),INTENT(IN)::coors
    INTEGER,DIMENSION(n_atoms_indices),INTENT(IN)::atoms_indices_serialized
    INTEGER,DIMENSION(n_atoms_series+1),INTENT(IN)::atoms_series_starts
    INTEGER,DIMENSION(n_frame_indices),INTENT(IN)::frame_indices
    DOUBLE PRECISION,DIMENSION(n_atoms_indices)::weights
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms_series,3)::center

    weights(:)=1.0d0

    center=CENTER_OF_MASS(coors, atoms_indices_serialized, atoms_series_starts, weights, frame_indices,&
                          n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frame_indices)

  END FUNCTION

  FUNCTION CENTER_OF_MASS(coors, atoms_indices_serialized, atoms_series_starts, weights, frame_indices,&
                          n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frame_indices) RESULT(center)

    IMPLICIT NONE

    INTEGER,INTENT(IN)::n_frames, n_atoms, n_atoms_indices, n_atoms_series, n_frame_indices
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms,3),INTENT(IN)::coors
    INTEGER,DIMENSION(n_atoms_indices),INTENT(IN)::atoms_indices_serialized
    INTEGER,DIMENSION(n_atoms_series+1),INTENT(IN)::atoms_series_starts
    DOUBLE PRECISION,DIMENSION(n_atoms_indices),INTENT(IN)::weights
    INTEGER,DIMENSION(n_frame_indices),INTENT(IN)::frame_indices
    DOUBLE PRECISION,DIMENSION(n_frames,n_atoms_series,3)::center

    INTEGER::ii,jj,ll,frame_index,atom_index
    DOUBLE PRECISION::total_weight

    center(:,:,:)=0.0d0

    DO ii=1,n_atoms_series
        total_weight=0.0d0
        DO jj=atoms_series_starts(ii)+1,atoms_series_starts(ii+1)
            total_weight = total_weight+weights(jj)
        END DO
        DO jj=atoms_series_starts(ii)+1,atoms_series_starts(ii+1)
            atom_index = atoms_indices_serialized(jj)+1
            DO ll=1,n_frame_indices
                frame_index=frame_indices(ll)+1
                center(ll,ii,:)=center(ll,ii,:)+weights(jj)*coors(frame_index,atom_index,:)/total_weight
            END DO
        END DO
    END DO

  END FUNCTION CENTER_OF_MASS

END MODULE MODULE_COM

max_stars_repo_name: LMMV/MolModMT

max_issues_repo_name: LMMV/MolModMT

max_forks_repo_name: LMMV/MolModMT


text:
      INTEGER FUNCTION RANMKT (STDDX,LAVSMT, REGDX,LAVRMT,IMPDX,LAVIMT)

C  RANMKT chooses a marker type at random from the three lists of
C  available marker types.

      INTEGER    STDDX,REGDX,IMPDX, WHLIST, RNDINT
      INTEGER    LAVSMT(*), LAVRMT(*), LAVIMT(*)

100   CONTINUE
      WHLIST = RNDINT (1,3)
      IF (WHLIST .EQ. 1) THEN
         IF (STDDX .LE. 0) GOTO 100
         RANMKT = LAVSMT(RNDINT(1,STDDX))
C don't allow marker type 1 (dot)
         IF (RANMKT .EQ. 1) GOTO 100
      ELSEIF (WHLIST .EQ. 2) THEN
         IF (REGDX .LE. 0) GOTO 100
         RANMKT = LAVRMT(RNDINT(1,REGDX))
      ELSE
         IF (IMPDX .LE. 0) GOTO 100
         RANMKT = LAVIMT(RNDINT(1,IMPDX))
      ENDIF

      END


text:
module quaff_physical_constants
    use quaff_acceleration_m, only: acceleration_t
    use quaff_energy_per_temperature_amount_m, only: energy_per_temperature_amount_t
    use quaff_pressure_m, only: pressure_t
    use quaff_conversion_factors_m, only: &
            GRAVITY_ => GRAVITY, &
            PASCALS_PER_ATMOSPHERE

    implicit none

    type(pressure_t), parameter :: ATMOSPHERIC_PRESSURE = &
            pressure_t(pascals = PASCALS_PER_ATMOSPHERE)
    type(acceleration_t), parameter :: GRAVITY = &
            acceleration_t(meters_per_square_second = GRAVITY_)
    type(energy_per_temperature_amount_t), parameter :: UNIVERSAL_GAS_CONSTANT = &
            energy_per_temperature_amount_t(joules_per_kelvin_mol = 8.31446261815324d0)
end module

max_stars_repo_name: everythingfunctional/quaff

max_issues_repo_name: everythingfunctional/quaff

max_forks_repo_name: everythingfunctional/quaff


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

module hecmw_precond_33
  use hecmw_util
  use hecmw_matrix_misc
  use hecmw_precond_BILU_33
  use hecmw_precond_DIAG_33
  use hecmw_precond_SSOR_33
  use hecmw_precond_ML_33
  use hecmw_precond_SAINV_33
  use hecmw_precond_RIF_33
  use hecmw_precond_nn
  use hecmw_solver_las_33
  implicit none

  private

  public :: hecmw_precond_33_setup
  public :: hecmw_precond_33_clear
  public :: hecmw_precond_33_apply

contains

  subroutine hecmw_precond_33_setup(hecMAT, hecMESH, sym)
    implicit none
    type (hecmwST_matrix), intent(inout) :: hecMAT
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in) :: sym

    select case(hecmw_mat_get_precond( hecMAT ))
      case(1,2)
        call hecmw_precond_SSOR_33_setup(hecMAT)
      case(3)
        call hecmw_precond_DIAG_33_setup(hecMAT)
      case(5)
        call hecmw_precond_ML_33_setup(hecMAT, hecMESH, sym)
      case(10,11,12)
        call hecmw_precond_BILU_33_setup(hecMAT)
      case(20)
        call hecmw_precond_33_SAINV_setup(hecMAT)
      case(21)
        call hecmw_precond_RIF_33_setup(hecMAT)
      case default
        call hecmw_precond_nn_setup(hecMAT, hecMESH, sym)
    end select

  end subroutine hecmw_precond_33_setup

  subroutine hecmw_precond_33_clear(hecMAT)
    implicit none
    type (hecmwST_matrix), intent(inout) :: hecMAT

    select case(hecmw_mat_get_precond( hecMAT ))
      case(1,2)
        call hecmw_precond_SSOR_33_clear(hecMAT)
      case(3)
        call hecmw_precond_DIAG_33_clear()
      case(5)
        call hecmw_precond_ML_33_clear()
      case(10:12)
        call hecmw_precond_BILU_33_clear()
      case(20)
        call hecmw_precond_33_SAINV_clear()
      case(21)
        call hecmw_precond_RIF_33_clear()
      case default
    end select

  end subroutine hecmw_precond_33_clear

  subroutine hecmw_precond_33_apply(hecMESH, hecMAT, R, Z, ZP, time_precond, COMMtime)
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix), intent(in)     :: hecMAT
    real(kind=kreal), intent(in) :: R(:)
    real(kind=kreal), intent(out) :: Z(:), ZP(:)
    real(kind=kreal), intent(inout) :: time_precond
    real(kind=kreal), intent(inout) :: COMMtime
    integer(kind=kint ) :: i, iterPRE, iterPREmax
    real(kind=kreal) :: START_TIME, END_TIME

    iterPREmax = hecmw_mat_get_iterpremax( hecMAT )
    do iterPRE= 1, iterPREmax
      START_TIME = hecmw_Wtime()
      select case(hecmw_mat_get_precond( hecMAT ))
        case(1,2)
          call hecmw_precond_SSOR_33_apply(ZP)
        case(3)
          call hecmw_precond_DIAG_33_apply(ZP)
        case(5)
          call hecmw_precond_ML_33_apply(ZP)
        case(10:12)
          call hecmw_precond_BILU_33_apply(ZP)
        case(20)
          call hecmw_precond_33_SAINV_apply(R,ZP)
        case(21)
          call hecmw_precond_RIF_33_apply(ZP)
        case default
      end select
      END_TIME = hecmw_Wtime()
      time_precond = time_precond + END_TIME - START_TIME

      !C-- additive Schwartz
      do i= 1, hecMAT%N * hecMAT%NDOF
        Z(i)= Z(i) + ZP(i)
      enddo
      if (iterPRE.eq.iterPREmax) exit

      !C--    {ZP} = {R} - [A] {Z}
      call hecmw_matresid_33 (hecMESH, hecMAT, Z, R, ZP, COMMtime)
    enddo
  end subroutine hecmw_precond_33_apply
end module hecmw_precond_33

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
*DECK CNBFA
      SUBROUTINE CNBFA (ABE, LDA, N, ML, MU, IPVT, INFO)
C***BEGIN PROLOGUE  CNBFA
C***PURPOSE  Factor a band matrix by elimination.
C***LIBRARY   SLATEC
C***CATEGORY  D2C2
C***TYPE      COMPLEX (SNBFA-S, DNBFA-D, CNBFA-C)
C***KEYWORDS  BANDED, LINEAR EQUATIONS, MATRIX FACTORIZATION,
C             NONSYMMETRIC
C***AUTHOR  Voorhees, E. A., (LANL)
C***DESCRIPTION
C
C    CNBFA factors a complex band matrix by elimination.
C
C     CNBFA is usually called by CNBCO, but it can be called
C     directly with a saving in time if RCOND is not needed.
C
C     On Entry
C
C        ABE     COMPLEX(LDA, NC)
C                contains the matrix in band storage.  The rows
C                of the original matrix are stored in the rows
C                of ABE and the diagonals of the original matrix
C                are stored in columns 1 through ML+MU+1 of ABE.
C                NC must be .GE. 2*ML+MU+1 .
C                See the comments below for details.
C
C        LDA     INTEGER
C                the leading dimension of the array ABE.
C                LDA must be .GE. N .
C
C        N       INTEGER
C                the order of the original matrix.
C
C        ML      INTEGER
C                number of diagonals below the main diagonal.
C                0 .LE. ML .LT. N .
C
C        MU      INTEGER
C                number of diagonals above the main diagonal.
C                0 .LE. MU .LT. N .
C                More efficient if ML .LE. MU .
C
C     On Return
C
C        ABE     an upper triangular matrix in band storage
C                and the multipliers which were used to obtain it.
C                the factorization can be written  A = L*U  where
C                L is a product of permutation and unit lower
C                triangular matrices and  U  is upper triangular.
C
C        IPVT    INTEGER(N)
C                an integer vector of pivot indices.
C
C        INFO    INTEGER
C                =0  normal value
C                =K  if  U(K,K) .EQ. 0.0 .  This is not an error
C                condition for this subroutine, but it does
C                indicate that CNBSL will divide by zero if
C                called.  Use RCOND in CNBCO for a reliable
C                indication of singularity.
C
C     Band Storage
C
C           If  A  is a band matrix, the following program segment
C           will set up the input.
C
C                   ML = (band width below the diagonal)
C                   MU = (band width above the diagonal)
C                   DO 20 I = 1, N
C                      J1 = MAX(1, I-ML)
C                      J2 = MIN(N, I+MU)
C                      DO 10 J = J1, J2
C                         K = J - I + ML + 1
C                         ABE(I,K) = A(I,J)
C                10    CONTINUE
C                20 CONTINUE
C
C           This uses columns  1  through  ML+MU+1  of ABE .
C           Furthermore,  ML  additional columns are needed in
C           ABE  starting with column  ML+MU+2  for elements
C           generated during the triangularization.  The total
C           number of columns needed in  ABE  is  2*ML+MU+1 .
C
C     Example:  If the original matrix is
C
C           11 12 13  0  0  0
C           21 22 23 24  0  0
C            0 32 33 34 35  0
C            0  0 43 44 45 46
C            0  0  0 54 55 56
C            0  0  0  0 65 66
C
C      then  N = 6, ML = 1, MU = 2, LDA .GE. 5  and ABE should contain
C
C            * 11 12 13  +     , * = not used
C           21 22 23 24  +     , + = used for pivoting
C           32 33 34 35  +
C           43 44 45 46  +
C           54 55 56  *  +
C           65 66  *  *  +
C
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED  CAXPY, CSCAL, CSWAP, ICAMAX
C***REVISION HISTORY  (YYMMDD)
C   800730  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CNBFA
      INTEGER LDA,N,ML,MU,IPVT(*),INFO
      COMPLEX ABE(LDA,*)
C
      INTEGER ML1,MB,M,N1,LDB,I,J,K,L,LM,LM1,LM2,MP,ICAMAX
      COMPLEX T
      COMPLEX ZDUM
      REAL CABS1
      CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM))
C
C***FIRST EXECUTABLE STATEMENT  CNBFA
      ML1=ML+1
      MB=ML+MU
      M=ML+MU+1
      N1=N-1
      LDB=LDA-1
      INFO=0
C
C     SET FILL-IN COLUMNS TO ZERO
C
      IF(N.LE.1)GO TO 50
      IF(ML.LE.0)GO TO 7
      DO 6 J=1,ML
        DO 5 I=1,N
          ABE(I,M+J)=(0.0E0,0.0E0)
    5   CONTINUE
    6 CONTINUE
    7 CONTINUE
C
C     GAUSSIAN ELIMINATION WITH PARTIAL ELIMINATION
C
      DO 40 K=1,N1
        LM=MIN(N-K,ML)
        LM1=LM+1
        LM2=ML1-LM
C
C     SEARCH FOR PIVOT INDEX
C
        L=-ICAMAX(LM1,ABE(LM+K,LM2),LDB)+LM1+K
        IPVT(K)=L
        MP=MIN(MB,N-K)
C
C     SWAP ROWS IF NECESSARY
C
        IF(L.NE.K)CALL CSWAP(MP+1,ABE(K,ML1),LDA,ABE(L,ML1+K-L),LDA)
C
C     SKIP COLUMN REDUCTION IF PIVOT IS ZERO
C
        IF(CABS1(ABE(K,ML1)).EQ.0.0E0) GO TO 20
C
C     COMPUTE MULTIPLIERS
C
        T=-(1.0E0,0.0E0)/ABE(K,ML1)
        CALL CSCAL(LM,T,ABE(LM+K,LM2),LDB)
C
C     ROW ELIMINATION WITH COLUMN INDEXING
C
        DO 10 J=1,MP
          CALL CAXPY (LM,ABE(K,ML1+J),ABE(LM+K,LM2),LDB,ABE(LM+K,LM2+J),
     1                LDB)
   10   CONTINUE
        GO TO 30
   20   CONTINUE
        INFO=K
   30   CONTINUE
   40 CONTINUE
   50 CONTINUE
      IPVT(N)=N
      IF(CABS1(ABE(N,ML1)).EQ.0.0E0) INFO=N
      RETURN
      END


text:
      SUBROUTINE CURVEO(XPT,YPT,LIMLO,LIMHI)
C
C          ------------------------------------------------
C          ROUTINE NO. (  92)   VERSION (A7.5)    11:FEB:85
C          ------------------------------------------------
C
C          THIS DRAWS AN OPEN CURVE THROUGH A SET OF POINTS.
C
C
C          THE ARGUMENTS ARE:
C
C          [XPT]    ARE THE X-COORDINATES AND
C          [YPT]    ARE THE Y-COORDINATES OF THE DATA POINTS,
C          <LIMLO>  IS THE LOWER LIMIT, AND
C          <LIMHI>  IS THE UPPER LIMIT OF THE POINTS IN THE ARRAYS.
C
C
      REAL XPT(LIMHI),YPT(LIMHI)
C
      COMMON /T0PPOS/ XPLOT0,YPLOT0
      COMMON /T0TRAC/ IPRINT
C
C
      IF (IPRINT.EQ.1) CALL G0MESG(35,0)
C
C          AUTOMATIC MAPPING WILL BE DONE IF NECESSARY.
C
      CALL G0AUTO(XPT,YPT,LIMLO,LIMHI,LIMLO,LIMHI,0)
C
      CALL G0CURV(XPT,YPT,LIMLO,LIMHI,0)
C
      XPLOT0= XPT(LIMHI)
      YPLOT0= YPT(LIMHI)
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
The Gifted Child Institute provides afterschool programs for gifted students in grade levels 2 through 6.  They also provide information and assistance to parents of gifted children.  Programs are currently offered at Robert E. Willett School Willett and Pioneer School Pioneer schools in Davis

Programs
  ProblemSolving
  Solve it. . .Engineer it!
  Solve it. . . Design it!


text:
      SUBROUTINE EQMCKA (IP,BGPDT,CSTM,EQEXIN,D,ISCALR)
C
C     ROUTINE FORMS D MATRIX (ACCTUALLY D TRANSPOSE)
C
      INTEGER         BGPDT,FILE,CSTM,EQEXIN,D,SYSBUF,IZ(1),MCB(7),
     1                NAME(2)
      REAL            TR(3,3),TI(3,3),DD(6,6),R(3)
      DIMENSION       TT(3,3)
C
      COMMON /SYSTEM/ SYSBUF
      COMMON /PACKX / IT1,IT2,II,JJ,INCR
      COMMON /ZZZZZZ/ Z(1)
C
      EQUIVALENCE     (IZ(1),Z(1))
C
      DATA    IZ2   , IZ3,IZ4,IZ5  / 2,3,4,5 /
      DATA    NAME  / 4HEQMC,4HKA  /
C
C     CONVERT  EXTERNAL IP TO INTERNAL IP
C
      IBUF = KORSZ(Z)-SYSBUF+1
      FILE = EQEXIN
      CALL GOPEN (EQEXIN,Z(IBUF),0)
      CALL READ  (*220,*10,EQEXIN,IZ(1),IBUF-1,0,IFLAG)
      GO TO 240
   10 CALL CLOSE (EQEXIN,1)
      DO 20 I = 1,IFLAG,2
      IF (IZ(I) .EQ.IP) GO TO 40
   20 CONTINUE
      CALL MESAGE (41,IP,NAME)
      IP = 0
      GO TO 50
   30 CALL MESAGE (41,IP,NAME)
C
C     SCALAR POINT
C
      GO TO 60
   40 IP = IZ(I+1)
C
C     FIND RZERO FOR  IP
C
   50 FILE = BGPDT
      R(1) = 0.0
      R(2) = 0.0
      R(3) = 0.0
      CALL GOPEN (BGPDT,Z(IBUF),0)
      IF (IP .EQ. 0) GO TO 70
      I= (IP-1)*4
      CALL FREAD (BGPDT,Z,-I,0)
      CALL FREAD (BGPDT,I, 1,0)
      IF (I .EQ. -1) GO TO 30
      CALL FREAD (BGPDT,R,3,0)
   60 CALL REWIND (BGPDT)
      CALL SKPREC (BGPDT,1)
C
C     SET UP TO WRITE D
C
   70 IBUF1 = IBUF-SYSBUF
      NZ = IBUF1-5
C
C     BRING IN CSTM
C
      FILE = CSTM
      CALL OPEN   (*90,CSTM,Z(IBUF1),0)
      CALL FWDREC (*220,CSTM)
      CALL READ   (*220,*80,CSTM,Z(IZ5),NZ,0,NCSTM)
      GO TO 240
   80 CALL CLOSE  (CSTM,1)
      CALL PRETRS (Z(IZ5),NCSTM)
   90 CALL GOPEN  (D,Z(IBUF1),1)
      CALL MAKMCB (MCB,D,6,2,1)
      ISCALR = 0
      II     = 1
      JJ     = 6
      IT1    = 1
      IT2    = 1
      INCR   = 1
C
C     EXAMINE BGPDT
C
  100 CALL READ (*220,*190,BGPDT,Z(1),4,0,IFLAG)
      IF (IZ(1) .EQ. -1) GO TO 170
C
C     COMPUTE  TR
C
      ISCALR  = 1
      TR(1,1) = 0.0
      TR(2,2) = 0.0
      TR(3,3) = 0.0
      TR(2,1) = Z(IZ4) -R(3)
      TR(1,2) =-TR(2,1)
      TR(3,1) = R(2)- Z(IZ3)
      TR(1,3) =-TR(3,1)
      TR(3,2) = Z(IZ2)-R(1)
      TR(2,3) =-TR(3,2)
      DO 110 I = 1,3
      DO 110 J = 1,3
      TI(I,J) = 0.0
      IF (I .EQ. J) TI(I,J) = 1.0
  110 CONTINUE
      IF (IZ(1) .EQ. 0) GO TO 130
      CALL TRANSS (IZ(1),TI)
      CALL GMMATS (TI,3,3,1,TR,3,3,0,TT)
      DO 120 I = 1,3
      DO 120 J = 1,3
  120 TR(I,J) = TT(I,J)
C
C     MOVE STUFF INTO  DD
C
  130 DO 150 I = 1,6
      DO 150 J = 1,3
      IF (I .GT. 3) GO TO 140
      DD(I  ,J  ) = TI(J,I)
      DD(I+3,J+3) = DD(I,J)
      GO TO 150
  140 DD(I,J) = TR(I-3,J)
      DD(J,I) = 0.0
  150 CONTINUE
      DO 160 I = 1,6
      CALL PACK (DD(1,I),D,MCB)
  160 CONTINUE
      GO TO 100
C
C     SCALAR POINT
C
  170 DO 180 I = 1,6
  180 DD(I,1) = 0.0
      CALL PACK (DD,D,MCB)
      GO TO 100
C
C     END BGPDT
C
  190 CALL CLOSE  (BGPDT,1)
      CALL CLOSE  (D,1)
      CALL WRTTRL (MCB)
      RETURN
C
C     ERROR MESAGES
C
  210 CALL MESAGE (IP1,FILE,NAME)
  220 IP1 = -2
      GO TO 210
  240 IP1 = -8
      GO TO 210
      END


text:
! Example of relativisation
!
! Stencil shapes are relative to the index on the left hand
! side. In this case the first index has a non-zero offset
! and so the right hand side is understood relative to this
program relativisation
implicit none

integer i, j, n
parameter (n = 3)
real a(0:n, 0:n)
real b(0:n, 0:n)

do i=1, n
   do j=1, n
      != stencil readOnce, backward(depth=2, dim=1, nonpointed)*pointed(dim=2) :: b
      a(i+1,j) = b(i,j) + b(i-1,j)
   end do
end do

end program

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


text:
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>c
c
      subroutine mpi_setup
c
c  Set up MPI stuff, including
c     - define the active set of processes
c     - set up new communicator
c
      include 'header.h'
c
      include 'mpi_stuff.h'
c
      integer no_nodes, color
c
c ... initialize MPI parameters
      call mpi_init(ierror)

      call mpi_comm_size(MPI_COMM_WORLD, no_nodes, ierror)
      call mpi_comm_rank(MPI_COMM_WORLD, myid, ierror)

      if (.not. convertdouble) then
         dp_type = MPI_DOUBLE_PRECISION
      else
         dp_type = MPI_REAL
      endif

c---------------------------------------------------------------------
c     let node 0 be the root for the group (there is only one)
c---------------------------------------------------------------------
      root = 0
c
      if (no_nodes .lt. num_procs) then
         if (myid .eq. root) write(*, 10) no_nodes, num_procs
   10    format(' Requested MPI processes ',i5,
     &          ' less than the compiled value ',i5)
         call mpi_abort(MPI_COMM_WORLD, 1, ierror)
         stop
      endif
c
      if (myid .ge. num_procs) then
         active = .false.
         color = 1
      else
         active = .true.
         color = 0
      end if

      call mpi_comm_split(MPI_COMM_WORLD,color,myid,comm_setup,ierror)
      if (.not. active) return

      call mpi_comm_rank(comm_setup, myid, ierror)
      if (no_nodes .ne. num_procs) then
         if (myid .eq. root) write(*, 20) no_nodes, num_procs
   20    format('Warning: Requested ',i5,' MPI processes, ',
     &          'but the compiled value is ',i5/
     &          'The compiled value is used for benchmarking')
      endif
c
      return
      end
c
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>c
c
      subroutine env_setup(tot_threads)
c
c  Set up from environment variables
c
c ... common variables
      include 'header.h'
c
      integer tot_threads
c
      include 'mpi_stuff.h'
c
c ... local variables
      integer ios, curr_threads, ip, mp, group, ip1, ip2
      integer entry_counts(num_procs)
      character envstr*80, line*132
c
c$    integer omp_get_max_threads
c$    external omp_get_max_threads
c
      if (myid .ne. root) goto 80
c
c ... test the OpenMP multi-threading environment
      mp = 0
c$    mp = omp_get_max_threads()
c
c ... master sets up parameters
      call getenv('OMP_NUM_THREADS', envstr)
      if (envstr .ne. ' ' .and. mp .gt. 0) then
         read(envstr,*,iostat=ios) num_threads
         if (ios.ne.0 .or. num_threads.lt.1) num_threads = 1
         if (mp .ne. num_threads) then
            write(*, 10) num_threads, mp
   10       format(' Warning: Requested ',i4,' threads per process,',
     &             ' but the active value is ',i4)
            num_threads = mp
         endif
      else
         num_threads = 1
      endif
c
      call getenv('NPB_MZ_BLOAD', envstr)
      if (envstr .ne. ' ') then
         if (envstr.eq.'on' .or. envstr.eq.'ON') then
            mz_bload = 1
         else if (envstr(1:1).eq.'t' .or. envstr(1:1).eq.'T') then
            mz_bload = 1
         else
            read(envstr,*,iostat=ios) mz_bload
            if (ios.ne.0) mz_bload = 0
         endif
      else
         mz_bload = 1
      endif
c
      call getenv('NPB_MAX_THREADS', envstr)
      max_threads = 0
      if (mz_bload.gt.0 .and. envstr.ne.' ') then
         read(envstr,*,iostat=ios) max_threads
         if (ios.ne.0 .or. max_threads.lt.0) max_threads = 0
         if (max_threads.gt.0 .and. max_threads.lt.num_threads) then
            write(*,20) max_threads, num_threads
   20       format(' Error: max_threads ',i5,
     &             ' is less than num_threads ',i5/
     &             ' Please redefine the value for NPB_MAX_THREADS',
     &             ' or OMP_NUM_THREADS')
            call mpi_abort(MPI_COMM_WORLD, 1, ierror)
            stop
         endif
      endif
c
      call getenv('NPB_VERBOSE', envstr)
      npb_verbose = 0
      if (envstr.ne.' ') then
         read(envstr,*,iostat=ios) npb_verbose
         if (ios.ne.0) npb_verbose = 0
      endif
c
      do ip = 1, num_procs
         proc_num_threads(ip) = num_threads
         proc_group(ip) = 0
      end do
c
      open(2, file='loadsp-mz.data', status='old', iostat=ios)
      if (ios.eq.0) then
         write(*,*) 'Reading load factors from loadsp-mz.data'

         if (mz_bload .ge. 1) then
            mz_bload = -mz_bload
         endif

         do ip = 1, num_procs
            entry_counts(ip) = 0
         end do

         do while (.true.)
   25       read(2,'(a)',end=40,err=40) line
            if (line.eq.' ' .or. line(1:1).eq.'#') goto 25

            call decode_line(line, ip1, ip2, curr_threads, group, ios)
            if (ios .ne. 0) goto 40

            if (mz_bload .lt. 0 .and. group .gt. 0) then
               mz_bload = -mz_bload
            endif

            if (curr_threads .lt. 1) curr_threads = 1
            if (mp .le. 0) curr_threads = 1
            if (ip1.lt.0) ip1 = 0
            if (ip2.ge.num_procs) ip2 = num_procs - 1

            do ip = ip1+1, ip2+1
               proc_num_threads(ip) = curr_threads
               proc_group(ip) = group
               entry_counts(ip) = entry_counts(ip) + 1
            end do
         end do
   40    close(2)

         do ip = 1, num_procs
            if (entry_counts(ip) .eq. 0) then
               write(*,*) '*** Error: Missing entry for proc ',ip-1
               call mpi_abort(MPI_COMM_WORLD, 1, ierror)
               stop
            else if (entry_counts(ip) .gt. 1) then
               write(*,*) '*** Warning: Multiple entries for proc ',
     &                    ip-1, ', only the last one used'
            endif
         end do

         ip1 = 1
      else
         write(*,*) 'Use the default load factors with threads'
         ip1 = 0
      endif

      if (ip1 .gt. 0 .or. npb_verbose .gt. 0) then
         ip1 = 0
         do ip = 1, num_procs
            if (ip .eq. 1 .or.
     &          proc_num_threads(ip) .ne. curr_threads .or.
     &          proc_group(ip) .ne. group) then

               ip2 = ip-2
               if (ip2 .gt. ip1+1) write(*,*) '    ...'
               if (ip2 .gt. ip1)
     &            write(*,30) ip2, curr_threads, group

               curr_threads = proc_num_threads(ip)
               group = proc_group(ip)

               ip1 = ip - 1
               write(*,30) ip1, curr_threads, group

            else if (ip .eq. num_procs) then
               ip2 = ip-1
               if (ip2 .gt. ip1+1) write(*,*) '    ...'
               write(*,30) ip2, curr_threads, group
            endif
         end do
   30    format('  proc',i6,'  num_threads =',i5,
     >          '  group =',i5)
      endif
c
c ... broadcast parameters to all processes
   80 call mpi_bcast(num_threads, 1, mpi_integer, root,
     &               comm_setup, ierror)
      call mpi_bcast(mz_bload, 1, mpi_integer, root,
     &               comm_setup, ierror)
      call mpi_bcast(max_threads, 1, mpi_integer, root,
     &               comm_setup, ierror)
      call mpi_bcast(proc_num_threads, num_procs, mpi_integer, root,
     &               comm_setup, ierror)
      call mpi_bcast(proc_group, num_procs, mpi_integer, root,
     &               comm_setup, ierror)
      call mpi_bcast(npb_verbose, 1, mpi_integer, root,
     &               comm_setup, ierror)
c
      tot_threads = 0
      do ip = 1, num_procs
         tot_threads = tot_threads + proc_num_threads(ip)
      end do
      if (myid .eq. root) then
         if (mp .gt. 0)
     &      write(*, 1004) tot_threads, dble(tot_threads)/num_procs
      endif
 1004 format(' Total number of threads: ', i6,
     &       '  (', f5.1, ' threads/process)')
c
      return
      end
c
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>c
c
      subroutine decode_line(line, ip1, ip2, curr_threads, group, ios)
      implicit none
c
c  decode a line from the load data file
c  format:  ip1[:ip2] curr_threads group
c
      character line*(*)
      integer ip1, ip2, curr_threads, group, ios
c
      integer is, n
c
      ios = -1
c
      n  = len(line)
      is = 1
      do while (is.le.n .and. line(is:is).ne.':')
         if (line(is:is).eq.'!') n = is
         is = is + 1
      end do
c
      if (is .gt. n) then
         read(line,*,err=90,end=90) ip1, curr_threads, group
         ip2 = ip1
      else if (is.eq.1 .or. is.eq.n) then
         go to 90
      else
         read(line(:is-1),*,err=90,end=90) ip1
         read(line(is+1:),*,err=90,end=90) ip2, curr_threads, group
      endif
c
      if (ip2 .lt. ip1) then
         is  = ip2
         ip2 = ip1
         ip1 = is
      endif
      ios = 0
c
   90 return
      end
c
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>c
c
      subroutine get_comm_index(zone, iproc, comm_index)
c
      include 'header.h'
c
      include 'mpi_stuff.h'
c
c  Calculate the communication index of a zone within a processor group
c
      integer zone, iproc, comm_index
c
c     local variables
      integer izone, jzone
c
      jzone  = (zone - 1)/x_zones + 1
      izone  = mod(zone - 1, x_zones) + 1
c
      comm_index = 0
      if (zone_proc_id(iz_west(zone)) .eq. iproc)
     $   comm_index = comm_index + y_size(jzone)
      if (zone_proc_id(iz_east(zone)) .eq. iproc)
     $   comm_index = comm_index + y_size(jzone)
      if (zone_proc_id(iz_south(zone)) .eq. iproc)
     $   comm_index = comm_index + x_size(izone)
      if (zone_proc_id(iz_north(zone)) .eq. iproc)
     $   comm_index = comm_index + x_size(izone)
c
      return
      end
c
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>c
c
      subroutine map_zones(num_zones, nx, ny, nz, tot_threads)
c
c  Perform zone-process mapping for load balance
c
      include 'header.h'
c
      integer num_zones, nx(*), ny(*), nz(*), tot_threads
c
      include 'mpi_stuff.h'
c
c     local variables
      integer z_order(max_zones)
      integer zone, iz, z2, mz, np, ip, zone_comm, comm_index
      integer n1, n2, work1, work2, np1, np2
      integer imx, imn, inc, icur_size
      double precision tot_size, cur_size, max_size, ave_size
      double precision zone_size(max_zones)
      double precision diff_ratio, tot_group_size
c
      integer group, ipg, tot_group_threads
      integer proc_group_flag(num_procs)
c
c ... sort the zones in decending order
      tot_size = 0.d0
      do iz = 1, num_zones
         zone_size(iz) = 1.0d0*nx(iz)*ny(iz)*nz(iz)
         z_order(iz) = iz
         tot_size = tot_size + zone_size(iz)
      end do
      do iz = 1, num_zones-1
         cur_size = zone_size(z_order(iz))
         mz = iz
         do z2 = iz+1, num_zones
            if (cur_size.lt.zone_size(z_order(z2))) then
               cur_size = zone_size(z_order(z2))
               mz = z2
            endif
         end do
         if (mz .ne. iz) then
            z2 = z_order(iz)
            z_order(iz) = z_order(mz)
            z_order(mz) = z2
         endif
      end do
c
      if (npb_verbose .gt. 1 .and. myid .eq. root) then
         write(*,10)
         do iz = 1, num_zones
            z2 = z_order(iz)
            write(*,15) iz,z2,nx(z2),ny(z2),nz(z2),zone_size(z2)
         end do
      endif
   10 format(/' Sorted zones:'/
     &       '  seq. zone    nx    ny    nz    size')
   15 format(i5,':',4(1x,i5),1x,f9.0)
c
c ... balance the load among processes
      do ip = 1, num_procs
         proc_zone_count(ip) = 0
         proc_zone_size(ip) = 0.d0
      end do

      if (abs(mz_bload).gt.1) goto 110
c
c ... try a simple block packing scheme
      n1 = 1
      n2 = num_procs
      work1 = mod(y_zones,num_procs)
      do while (n1 .le. n2)
         if (n1*n2 .eq. num_procs) then
            work2 = mod(x_zones,n1) + mod(y_zones,n2)
            if (work2 .le. work1) then
               work1 = work2
               np1 = n1
               np2 = n2
            endif
         endif
         n1 = n1 + 1
         n2 = num_procs / n1
      end do
c ... if can't find a good solution, fall back to next method
      if (work1 .ne. 0) goto 110
c
c ... amount of work for each block
      work1 = x_zones / np1
      work2 = y_zones / np2
c
      do iz = 1, num_zones
         n1 = mod(iz-1,x_zones) / work1
         n2 = (iz - 1) / x_zones / work2
         ip = n1 + n2*np1 + 1
c
c  ...   assign the zone to the current processor group
         zone_proc_id(iz) = ip - 1
         proc_zone_size(ip) = proc_zone_size(ip) + zone_size(iz)
         proc_zone_count(ip) = proc_zone_count(ip) + 1
      end do
c
      goto 300
c
c ... use a bin-packing scheme to balance the load among processes
  110 continue
      do iz = 1, num_zones
         zone_proc_id(iz) = -1
      end do

      iz = 1
      do while (iz .le. num_zones)
c
c  ...   the current most empty processor
         np = 1
         cur_size = proc_zone_size(1)
         do ip = 2, num_procs
            if (cur_size.gt.proc_zone_size(ip)) then
               np = ip
               cur_size = proc_zone_size(ip)
            endif
         end do
         ip = np - 1
c
c  ...   get a zone that has the largest communication index with
c        the current group and does not worsen the computation balance
         mz = z_order(iz)
         if (iz .lt. num_zones) then
            call get_comm_index(mz, ip, zone_comm)
            do z2 = iz+1, num_zones
               zone = z_order(z2)

               diff_ratio = (zone_size(z_order(iz)) -
     &                      zone_size(zone)) / zone_size(z_order(iz))
               if (diff_ratio .gt. 0.05D0) goto 120

               if (zone_proc_id(zone) .lt. 0) then
                  call get_comm_index(zone, ip, comm_index)
                  if (comm_index .gt. zone_comm) then
                     mz = zone
                     zone_comm = comm_index
                  endif
               endif
            end do
         endif
c
c  ...   assign the zone to the current processor group
  120    zone_proc_id(mz) = ip
         proc_zone_size(np) = proc_zone_size(np) + zone_size(mz)
         proc_zone_count(np) = proc_zone_count(np) + 1
c
c  ...   skip the previously assigned zones
         do while (iz.le.num_zones)
            if (zone_proc_id(z_order(iz)).lt.0) goto 130
            iz = iz + 1
         end do
  130    continue
      end do
c
c ... move threads around if needed
      mz = 1
      if (tot_threads.eq.num_procs .or. mz_bload.lt.1) mz = 0
c
      if (mz .ne. 0) then
c
         do ipg = 1, num_procs
            proc_group_flag(ipg) = 0
         end do
c
         ipg = 1
c
c ...    balance load within a processor group
  200    do while (ipg .le. num_procs)
            if (proc_group_flag(ipg) .eq. 0) goto 210
            ipg = ipg + 1
         end do
  210    if (ipg .gt. num_procs) goto 300
c
         group = proc_group(ipg)
         tot_group_size = 0.d0
         tot_group_threads = 0
         do ip = ipg, num_procs
            if (proc_group(ip) .eq. group) then
               proc_group_flag(ip) = 1
               tot_group_size = tot_group_size + proc_zone_size(ip)
               tot_group_threads = tot_group_threads +
     &                             proc_num_threads(ip)
            endif
         end do
c
         ave_size = tot_group_size / tot_group_threads
c
c  ...   distribute size evenly among threads
         icur_size = 0
         do ip = 1, num_procs
            if (proc_group(ip) .ne. group) goto 220
            proc_num_threads(ip) = proc_zone_size(ip) / ave_size
            if (proc_num_threads(ip) .lt. 1)
     &          proc_num_threads(ip) = 1
            if (max_threads .gt. 0 .and.
     &          proc_num_threads(ip) .gt. max_threads)
     &          proc_num_threads(ip) = max_threads
            icur_size = icur_size + proc_num_threads(ip)
  220       continue
         end do
         mz = tot_group_threads - icur_size
c
c  ...   take care of any remainers
         inc = 1
         if (mz .lt. 0) inc = -1
         do while (mz .ne. 0)
            max_size = 0.d0
            imx = 0
            do ip = 1, num_procs
               if (proc_group(ip) .ne. group) goto 230
               if (mz .gt. 0) then
                  cur_size = proc_zone_size(ip) / proc_num_threads(ip)
                  if (cur_size.gt.max_size .and. (max_threads.le.0
     &                .or. proc_num_threads(ip).lt.max_threads)) then
                     max_size = cur_size
                     imx = ip
                  endif
               else if (proc_num_threads(ip) .gt. 1) then
                  cur_size = proc_zone_size(ip) /
     &                       (proc_num_threads(ip)-1)
                  if (max_size.eq.0 .or. cur_size.lt.max_size) then
                     max_size = cur_size
                     imx = ip
                  endif
               endif
  230          continue
            end do
            proc_num_threads(imx) = proc_num_threads(imx) + inc
            mz = mz - inc
         end do
c
         goto 200
      endif
c
c ... print the mapping
  300 if (npb_verbose .gt. 0 .and. myid .eq. root) then
         write(*,20)
         do ip = 1, num_procs
            write(*,25) ip-1,proc_zone_count(ip),
     &            proc_zone_size(ip),proc_num_threads(ip),
     &            proc_zone_size(ip)/proc_num_threads(ip)
            do iz = 1, num_zones
               if (zone_proc_id(iz) .eq. ip-1) then
                  write(*,30) iz, zone_size(iz)
               endif
            end do
         end do
      endif
   20 format(/' Zone-process mapping:'/
     &       '  proc nzones  zone_size nthreads size_per_thread')
   25 format(i5,2x,i5,2x,f10.0,2x,i5,3x,f10.0)
   30 format(3x,'zone ',i5,2x,f9.0)
c
      if (myid .eq. root) then
         imx = 1
         max_size = proc_zone_size(1)/proc_num_threads(1)
         imn = imx
         ave_size = max_size
         do ip = 2, num_procs
            cur_size = proc_zone_size(ip)/proc_num_threads(ip)
            if (cur_size.gt.max_size) then
               imx = ip
               max_size = cur_size
            endif
            if (cur_size.lt.ave_size) then
               imn = ip
               ave_size = cur_size
            endif
         end do

         if (npb_verbose .gt. 0) then
            write(*,*)
            write(*,35) 'Max', imx-1, proc_zone_count(imx),
     &                  proc_zone_size(imx),proc_num_threads(imx)
            write(*,35) 'Min', imn-1, proc_zone_count(imn),
     &                  proc_zone_size(imn),proc_num_threads(imn)
         endif
   35    format(1x,a,': proc=',i5,' nzones=',i5,' size=',f10.0,
     &          ' nthreads=',i5)

         write(*,40) tot_size / max_size
      endif
   40 format(/' Calculated speedup = ',f9.2/)
c
c ... reorganize list of zones for this process
      zone = 0
      do iz = 1, num_zones
         if (zone_proc_id(iz) .eq. myid) then
            zone = zone + 1
            proc_zone_id(zone) = iz
         endif
      end do
      proc_num_zones = zone
      if (zone .ne. proc_zone_count(myid+1)) then
         write(*,*) 'Warning: ',myid, ': mis-matched zone counts -',
     &              zone, proc_zone_count(myid+1)
      endif
c
c ... set number of threads for this process
      group = proc_group(myid+1)
      np = 0
      do ip = 1, num_procs
         if (proc_group(ip) .eq. group) then
            proc_group(ip) = np
            np = np + 1
            proc_num_threads(np) = proc_num_threads(ip)
         endif
      end do
      ipg = proc_group(myid+1)
      if (npb_verbose.gt.1) then
         write(*,50) myid, group, np, ipg, proc_num_threads(ipg+1)
      endif
   50 format(' myid',i5,' group',i5,' group_size',i5,
     &       ' group_pid',i5,' threads',i4)
c$    call omp_set_num_threads(proc_num_threads(ipg+1))
c
c ... pin-to-node within one process group
c      call smp_pinit_thread(np, ipg, proc_num_threads)
c
      return
      end

max_stars_repo_name: robert-mijakovic/readex-ptf

max_issues_repo_name: robert-mijakovic/readex-ptf

max_forks_repo_name: robert-mijakovic/readex-ptf


text:
if (myid==0) write(*,*) "nn_PML=", nn_PML
if (nn_PML > 0) then
    if (myid==0) write(*,*) "reallocating p, w, dg"
    deallocate(p)
    deallocate(w)
    deallocate(dg)
    allocate(p(ndf,nn+nn_PML),stat=error_id)
    allocate(w(ndf,nn+nn_PML),stat=error_id)
    allocate(dg(ndf,nn+nn_PML),stat=error_id)
else
    if (myid==0) write(*,*) "no reallocation of p, w, dg necessary"
endif

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
The Seminars Project is an archive of recorded lectures and seminars from the various University of California UC campuses collected and made available to the public by UCTV.  The project is the brain child of UC Davis entomology professor James R. Carey, who encouraged the recording, broadcast, and archival of many of the hundreds of seminars that take place on a weekly basis throughout the UC System.

The Enterprise http://wwwnew.davisenterprise.com/localnews/ucd/ucseminarsonlinefromecologytoevolutiontoelectronicteaching/ reported that an initial pilot program local to UC Davis made 49 seminars available to the public and met with a great deal of approval.  As of May 2011, seminars are available from the UC campuses at Berkeley, Davis, Los Angeles, San Diego, San Francisco, Santa Barbara, and the Washington, D.C. extension.

The Projects seminars are available at http://seminars.uctv.tv/ where they can be browsed by campus, by conference, by host (department, group, or organization, apparently), by speaker, or by subject matter.  Broad subjects include Arts & Humanities, Engineering, Life Sciences and Health, Physical Sciences & Mathematics, and Social & Behavioral Sciences, with subcategories in each.


text:
C     The two arrays are not overlapped

      PROGRAM ALIAS
      REAL A(10000)
      M = 50
      CALL FOO(A,M)
      END

      SUBROUTINE FOO(X,N)
      REAL X(100,100)
      L = 20
      DO I=1,L
         DO J =1,L
            CALL TOTO(X(I,J),X(2*I,2*J),I,J,L)
         ENDDO
      ENDDO
      END

      SUBROUTINE TOTO(V1,V2,II,JJ,L)
      REAL V1(L),V2(L)
      DO I=1,L
         V1(I)=I
      ENDDO
      END


subroutine locate(XX,N,X,J)
  implicit none
  !
  ! Integers
  integer*4 J,JL,JM,JU,N
  !
  ! Floats
  real*8 X,XX(N)
  !
  ! Code
  JL=0
  JU=N+1
  do while (JU-JL.gt.1)
     JM=(JU+JL)/2
     if ((XX(N).gt.XX(1)).eqv.(X.gt.XX(JM))) then
        JL=JM
     else
        JU=JM
     endif
  end do
  J=JL
  return
end subroutine locate