1*
2* $Id$
3*
4
5      subroutine kerker_G_init()
6      implicit none
7
8#include "bafdecls.fh"
9#include "errquit.fh"
10
11*     **** local variables ****
12      integer npack0,nfft3d,G(3)
13      integer i,ccode
14      real*8  gg,gg0,g0
15      integer tmp1(2)
16
17*     **** external functions ****
18      real*8   control_kerker_g0
19      external control_kerker_g0
20      integer  G_indx,c_G_indx,control_code
21      external G_indx,c_G_indx,control_code
22
23*     **** common block used for kerker_G.F ****
24c     real*8 tg(nfft3d)
25      logical isband,dokerker
26      integer tg_indx,tg_hndl
27      common / kerker_g_block / tg_indx,tg_hndl,isband,dokerker
28
29
30      g0 = control_kerker_g0()
31
32!$OMP MASTER
33      dokerker =  (g0.gt.0.0d0)
34!$OMP END MASTER
35!$OMP BARRIER
36
37      if (dokerker) then
38
39      ccode = control_code()
40!$OMP MASTER
41      isband =((ccode.eq.5).or.(ccode.eq.13).or.(ccode.eq.14))
42!$OMP END MASTER
43!$OMP BARRIER
44
45      if (isband) then
46         call C3dB_nfft3d(1,nfft3d)
47         call Cram_npack(0,npack0)
48         G(1)= c_G_indx(1)
49         G(2)= c_G_indx(2)
50         G(3)= c_G_indx(3)
51      else
52         call D3dB_nfft3d(1,nfft3d)
53         call Pack_npack(0,npack0)
54         G(1)= G_indx(1)
55         G(2)= G_indx(2)
56         G(3)= G_indx(3)
57      end if
58
59      if (.not.BA_alloc_get(mt_dbl,npack0,
60     >                     'tg',tg_hndl,tg_indx))
61     >   call errquit('kerker_G_init:out of heap memory',0,MA_ERR)
62
63      if (.not.BA_push_get(mt_dbl,nfft3d,'tmp1',tmp1(2),tmp1(1)))
64     >   call errquit('kerker_G_init:out of stack memory',0,MA_ERR)
65
66
67      gg0 = g0*g0
68!$OMP DO
69      do i = 1,nfft3d
70         gg  = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1)
71     >         + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1)
72     >         + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1))
73         dbl_mb(tmp1(1)+i-1) = gg/(gg+gg0)
74      end do
75!$OMP END DO
76
77      if (isband) then
78         call Cram_r_pack(0,dbl_mb(tmp1(1)))
79         call Cram_r_Copy(0,dbl_mb(tmp1(1)),dbl_mb(tg_indx))
80      else
81         call Pack_t_pack(0,dbl_mb(tmp1(1)))
82         call Pack_t_Copy(0,dbl_mb(tmp1(1)),dbl_mb(tg_indx))
83      end if
84      if (.not.BA_pop_stack(tmp1(2)))
85     >   call errquit('kerker_G_init:popping stack memory',0,MA_ERR)
86
87      end if
88      return
89      end
90
91      subroutine kerker_G_end()
92      implicit none
93
94#include "bafdecls.fh"
95#include "errquit.fh"
96
97*     **** common block used for kerker_G.F ****
98      logical isband,dokerker
99      integer tg_indx,tg_hndl
100      common / kerker_G_block / tg_indx,tg_hndl,isband,dokerker
101
102      if (dokerker) then
103         if (.not.BA_free_heap(tg_hndl))
104     >   call errquit('error freeing heap',0, MA_ERR)
105      end if
106      return
107      end
108
109
110      subroutine kerker_G(V)
111      implicit none
112      real*8     V(*)
113
114#include "bafdecls.fh"
115#include "errquit.fh"
116
117*     **** common block used for kerker_G.F ****
118      logical isband,dokerker
119      integer tg_indx,tg_hndl
120      common / kerker_g_block / tg_indx,tg_hndl,isband,dokerker
121
122*     **** local variables ****
123      integer nfft3d,n2ft3d,tmp1(2)
124      integer nx,ny,nz
125      real*8  scal1
126
127      if (dokerker) then
128
129      if (isband) then
130         call C3dB_nx(1,nx)
131         call C3dB_ny(1,ny)
132         call C3dB_nz(1,nz)
133         call C3dB_nfft3d(1,nfft3d)
134         n2ft3d = nfft3d
135      else
136         call D3dB_nx(1,nx)
137         call D3dB_ny(1,ny)
138         call D3dB_nz(1,nz)
139         call D3dB_nfft3d(1,nfft3d)
140         n2ft3d = 2*nfft3d
141      end if
142      scal1 = 1.0d0/dble(nx*ny*nz)
143
144
145      if (.not.BA_push_get(mt_dcpl,nfft3d,'tmp1',tmp1(2),tmp1(1)))
146     >   call errquit('kerker_G: out of stack memory',0, MA_ERR)
147      if (isband) then
148         call C3dB_rc_SMul(1,scal1,V,dcpl_mb(tmp1(1)))
149         call C3dB_rc_fft3f(1,dcpl_mb(tmp1(1)))
150         call Cram_c_pack(0,dcpl_mb(tmp1(1)))
151         call Cram_rc_Mul2(0,dbl_mb(tg_indx),dcpl_mb(tmp1(1)))
152         call Cram_c_unpack(0,dcpl_mb(tmp1(1)))
153         call C3dB_cr_fft3b(1,dcpl_mb(tmp1(1)))
154         call C3dB_cr_real(1,dcpl_mb(tmp1(1)),V)
155      else
156         call D3dB_r_SMul(1,scal1,V,dcpl_mb(tmp1(1)))
157         call D3dB_r_Zero_Ends(1,dcpl_mb(tmp1(1)))
158         call D3dB_rc_fft3f(1,dcpl_mb(tmp1(1)))
159         call Pack_c_pack(0,dcpl_mb(tmp1(1)))
160         call Pack_tc_Mul2(0,dbl_mb(tg_indx),dcpl_mb(tmp1(1)))
161         call Pack_c_unpack(0,dcpl_mb(tmp1(1)))
162         call D3dB_cr_fft3b(1,dcpl_mb(tmp1(1)))
163         call D3dB_r_Zero_Ends(1,dcpl_mb(tmp1(1)))
164         !call dcopy(n2ft3d,dcpl_mb(tmp1(1)),1,V,1)
165         call Parallel_shared_vector_copy(.true.,n2ft3d,
166     >                dcpl_mb(tmp1(1)),V)
167      end if
168      if (.not.BA_pop_stack(tmp1(2)))
169     >   call errquit('kerker_G: popping stack memory',0,MA_ERR)
170
171
172      endif
173      return
174      end
175
176