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