Branch data Line data Source code
1 : : !!****m* ABINIT/m_opernla_ylm_mv
2 : : !! NAME
3 : : !! m_opernla_ylm_mv
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2021-2026 ABINIT group (LB,MT)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_opernla_ylm_mv
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_errors
26 : : use m_xmpi
27 : : #if defined HAVE_OPENMP
28 : : use OMP_LIB
29 : : #endif
30 : :
31 : : use defs_abitypes, only : MPI_type
32 : : use m_time, only : timab
33 : :
34 : : implicit none
35 : :
36 : : private
37 : : !!***
38 : :
39 : : public :: opernla_ylm_mv
40 : : integer,public,save :: opernla_mv_counter = -1
41 : : integer,public,save :: opernla_mv_dgemv_counter = -1
42 : : !!***
43 : : !!***
44 : :
45 : : contains
46 : : !!***
47 : :
48 : : !!****f* ABINIT/opernla_ylm_mv
49 : : !! NAME
50 : : !! opernla_ylm_mv
51 : : !!
52 : : !! FUNCTION
53 : : !! "matrix-vector" alternative implementation of "opernla_ylm".
54 : : !!
55 : : !! For a given wave-function |c>, get all projected scalars
56 : : !! <p_lmn|c> where |p_lmn> are non-local projectors
57 : : !! With:
58 : : !! <p_lmn|c>=4pi/sqrt(vol) (i)^l Sum_g[c(g).f_nl(g).Y_lm(g).exp(2pi.i.g.R)]
59 : : !!
60 : : !! Here this is done in 3 steps:
61 : : !! (1) compute for every g : scal(g) = c(g).exp(2pi.i.g.R)
62 : : !! (2) compute for every lmn : scal(lmn) = Sum_g[scal(g).f_nl(g).Y_lm(g)]
63 : : !! (3) compute for every lmn : <p_lmn|c> = 4pi/sqrt(vol).(i)^l.scal(lmn)
64 : : !!
65 : : !! Step (2) is a real-matrix/complex-vector multiplication, here two options are possible:
66 : : !! - case nloalg(1)=2 : compute the real and imaginary parts separately using two calls of DGMEV
67 : : !! - case nloalg(1)=3 : in order to read the matrix only once, we compute both real and imaginary parts at the same time "by hand"
68 : : !!
69 : : !! Depending on the achitecture and the available blas library, one option could be more interesting than an other...
70 : : !!
71 : : !! INPUTS
72 : : !! choice=chooses possible output:
73 : : !! if choice>=0: compute projected scalars
74 : : !! if choice<0: same as choice>0 but use already computed projected scalars
75 : : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
76 : : !! 2 if <p_lmn|c> scalars are complex
77 : : !! dimffnl=second dimension of ffnl
78 : : !! ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
79 : : !! ia3=gives the number of the first atom in the subset presently treated
80 : : !! indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
81 : : !! istwf_k=option parameter that describes the storage of wfs
82 : : !! matblk=dimension of the array ph3d
83 : : !! mpi_enreg=information about MPI parallelization
84 : : !! nincat=number of atoms in the subset here treated
85 : : !! nlmn=number of (l,m,n) numbers for current type of atom
86 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
87 : : !! npw=number of plane waves in reciprocal space
88 : : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
89 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
90 : : !! ucvol=unit cell volume (bohr^3)
91 : : !! vect(2,npw*my_nspinor)=starting vector in reciprocal space
92 : : !!
93 : : !! OUTPUT
94 : : !! gx(cplex,nlmn,nincat,nspinor)= projected scalars
95 : : !!
96 : : !! SIDE EFFECTS
97 : : !!
98 : : !! NOTES
99 : : !! 1-Not available yet for openMP
100 : : !! 2-Operate for one type of atom, and within this given type of atom,
101 : : !! for a subset of at most nincat atoms.
102 : : !! 3-projector derivatives (abs(choice)>1) are not implemented yet
103 : : !!
104 : : !! SOURCE
105 : :
106 : 0 : subroutine opernla_ylm_mv(choice,cplex,dimffnl,ffnl,gx,&
107 : 0 : & ia3,indlmn,istwf_k,matblk,mpi_enreg,nincat,nlmn,&
108 [ # # ]: 0 : & nloalg,npw,nspinor,ph3d,ucvol,vect)
109 : :
110 : : !Arguments ------------------------------------
111 : : !scalars
112 : : integer,intent(in) :: choice,cplex,dimffnl,ia3,istwf_k,matblk
113 : : integer,intent(in) :: nincat,nlmn,npw,nspinor
114 : : real(dp),intent(in) :: ucvol
115 : : type(MPI_type),intent(in) :: mpi_enreg
116 : : !arrays
117 : : integer,intent(in) :: indlmn(6,nlmn),nloalg(3)
118 : : real(dp),intent(in),target :: ffnl(npw,dimffnl,nlmn)
119 : : real(dp),intent(in) :: ph3d(2,npw,matblk)
120 : : real(dp),intent(in) :: vect(:,:)
121 : : real(dp),intent(out) :: gx(cplex,nlmn,nincat,nspinor)
122 : :
123 : : !Local variables-------------------------------
124 : : !scalars
125 : : logical :: use_dgemv
126 : : integer :: ia,iaph3d,ierr,il,ilmn,ipw,ipw0,ipwshft,ispinor,jpw
127 : : real(dp) :: wt
128 : : !arrays
129 : : real(dp) :: buffer_r,buffer_i,tsec(2)
130 : 0 : real(dp),pointer :: ffnl_loc(:,:)
131 : 0 : real(dp),allocatable :: scali(:),scalr(:)
132 : 0 : real(dp),allocatable :: scalr_lmn(:),scali_lmn(:)
133 : : complex(dp) :: ctmp,cil(4)
134 : : ! *************************************************************************
135 : :
136 : 0 : if (choice==-1) return
137 : :
138 : : !Some checks
139 [ # # ]: 0 : if (abs(choice)>1) then
140 : 0 : ABI_ERROR('Only abs(choice)<=1 is available for now.')
141 : : end if
142 [ # # ]: 0 : if (nloalg(1)<2.or.nloalg(1)>10) then
143 : 0 : ABI_ERROR('nloalg(1) should be between 2 and 10.')
144 : : end if
145 : : ! nthreads=1
146 : : !#if defined HAVE_OPENMP
147 : : ! nthreads=OMP_GET_NUM_THREADS()
148 : : !#endif
149 : : ! if (nthreads>1) then
150 : : ! ABI_ERROR('Only nthreads=1 is available for now.')
151 : : ! end if
152 : :
153 [ # # ]: 0 : use_dgemv = nloalg(1)==2.or.nloalg(1)==5.or.nloalg(1)==7
154 : : if (choice>=0.or.abs(choice)>1) then
155 [ # # ]: 0 : if (use_dgemv) then
156 [ # # ]: 0 : if(opernla_mv_dgemv_counter>=0) opernla_mv_dgemv_counter = opernla_mv_dgemv_counter + 1
157 : : else
158 [ # # ]: 0 : if(opernla_mv_counter>=0) opernla_mv_counter = opernla_mv_counter + 1
159 : : end if
160 : : end if
161 : :
162 : : !Useful variables
163 [ # # ]: 0 : wt=four_pi/sqrt(ucvol);if (cplex==1) wt=2.d0*wt
164 [ # # # # ]: 0 : ipw0=1;if (istwf_k==2.and.mpi_enreg%me_g0_fft==1) ipw0=2
165 : :
166 : : !Allocate work space
167 [ # # # # : 0 : ABI_MALLOC(scalr,(npw))
# # # # #
# # # ]
168 [ # # # # : 0 : ABI_MALLOC(scali,(npw))
# # # # ]
169 [ # # # # : 0 : ABI_MALLOC(scalr_lmn,(nlmn))
# # # # #
# # # ]
170 [ # # # # : 0 : ABI_MALLOC(scali_lmn,(nlmn))
# # # # ]
171 : :
172 : 0 : ffnl_loc => ffnl(:,1,:)
173 : :
174 : : ! i^l
175 : 0 : cil(1) = ( 1.0_DP, 0.0_DP) * wt
176 : 0 : cil(2) = ( 0.0_DP, 1.0_DP) * wt
177 : 0 : cil(3) = (-1.0_DP, 0.0_DP) * wt
178 : 0 : cil(4) = ( 0.0_DP,-1.0_DP) * wt
179 : :
180 : : ! FIXME OpenMP parallelism is still flawed here so forcing NUM_THREADS=1 here
181 : : ! Bug reproducible with v9[71] and bounds checks enabled
182 : : !$OMP PARALLEL PRIVATE(il,ilmn,ipw,jpw), &
183 : : !$OMP PRIVATE(ispinor,ipwshft,ia,iaph3d) &
184 : : !$OMP NUM_THREADS(1)
185 : :
186 : : !Loop on spinorial components
187 [ # # ]: 0 : do ispinor =1,nspinor
188 : 0 : ipwshft=(ispinor-1)*npw
189 : :
190 : : ! Loop on atoms (blocking)
191 [ # # ]: 0 : do ia=1,nincat
192 [ # # ]: 0 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
193 : : ! Step (1) : Compute scal(g) = c(g).exp(2pi.i.g.R)
194 : : !$OMP DO
195 [ # # ]: 0 : do ipw=ipw0,npw
196 : 0 : jpw=ipw+ipwshft
197 : 0 : scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d))
198 : 0 : scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d))
199 : : end do
200 : : !$OMP END DO
201 : :
202 : : !$OMP SINGLE
203 [ # # ]: 0 : if (ipw0==2) then
204 : 0 : scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d)
205 : 0 : scali(1)=half*vect(1,1+ipwshft)*ph3d(2,1,iaph3d)
206 : : end if
207 : : !$OMP END SINGLE
208 : :
209 : : ! --------------------------------------------------------------------
210 : : ! ALL CHOICES:
211 : : ! Accumulate Gx
212 : : ! --------------------------------------------------------------------
213 : :
214 [ # # ]: 0 : if (choice>=0) then
215 : :
216 : : ! Step (2) : Compute scal(lmn) = Sum_g[scal(g).f_nl(g).Y_lm(g)]
217 [ # # ]: 0 : if (use_dgemv) then
218 [ # # ]: 0 : call DGEMV('T',npw,nlmn,1.0_DP,ffnl_loc,npw,scalr,1,0.0_DP,scalr_lmn,1)
219 [ # # ]: 0 : call DGEMV('T',npw,nlmn,1.0_DP,ffnl_loc,npw,scali,1,0.0_DP,scali_lmn,1)
220 : : else
221 [ # # ]: 0 : do ilmn=1,nlmn
222 : : ! do ipw=1,npw
223 : : ! scalr_lmn(ilmn) = scalr_lmn(ilmn) + scalr(ipw) * ffnl_loc(ipw,ilmn)
224 : : ! scali_lmn(ilmn) = scali_lmn(ilmn) + scali(ipw) * ffnl_loc(ipw,ilmn)
225 : : ! end do
226 : : !$OMP SINGLE
227 : : buffer_r = 0.0_DP
228 : : buffer_i = 0.0_DP
229 : : !$OMP END SINGLE
230 : : !$OMP DO REDUCTION(+:buffer_r,buffer_i)
231 [ # # ]: 0 : do ipw=1,npw
232 : 0 : buffer_r = buffer_r + scalr(ipw) * ffnl_loc(ipw,ilmn)
233 : 0 : buffer_i = buffer_i + scali(ipw) * ffnl_loc(ipw,ilmn)
234 : : end do
235 : : !$OMP END DO
236 : : !$OMP SINGLE
237 : 0 : scalr_lmn(ilmn) = buffer_r
238 : 0 : scali_lmn(ilmn) = buffer_i
239 : : !$OMP END SINGLE
240 : : end do
241 : : end if
242 : : ! Step (3) : Compute gx(lmn) = 4pi/sqrt(vol) (i)^l scal(lmn)
243 : : !$OMP SINGLE
244 [ # # ]: 0 : if (cplex==2) then
245 [ # # ]: 0 : do ilmn=1,nlmn
246 : 0 : il=mod(indlmn(1,ilmn),4)+1
247 : 0 : ctmp = cil(il) * cmplx(scalr_lmn(ilmn),scali_lmn(ilmn),kind=DP)
248 : 0 : gx(1,ilmn,ia,ispinor) = real(ctmp)
249 : 0 : gx(2,ilmn,ia,ispinor) = aimag(ctmp)
250 : : end do
251 : : else
252 [ # # ]: 0 : do ilmn=1,nlmn
253 : 0 : il=mod(indlmn(1,ilmn),4)+1
254 : 0 : ctmp = cil(il) * cmplx(scalr_lmn(ilmn),scali_lmn(ilmn),kind=DP)
255 : 0 : gx(1,ilmn,ia,ispinor) = real(ctmp)
256 : : end do
257 : : end if
258 : : !$OMP END SINGLE
259 : :
260 : : end if
261 : :
262 : : end do ! End loop on atoms
263 : :
264 : : end do ! End loop on spinorial components
265 : : !$OMP END PARALLEL
266 : :
267 : : !Deallocate temporary space
268 [ # # ]: 0 : ABI_FREE(scalr)
269 [ # # ]: 0 : ABI_FREE(scali)
270 [ # # ]: 0 : ABI_FREE(scalr_lmn)
271 [ # # ]: 0 : ABI_FREE(scali_lmn)
272 : :
273 : : !Has to reduce arrays in case of FFT parallelization
274 [ # # ]: 0 : if (mpi_enreg%nproc_fft>1) then
275 : 0 : call timab(48,1,tsec)
276 [ # # ]: 0 : if (choice>=0) then
277 : 0 : call xmpi_sum(gx,mpi_enreg%comm_fft,ierr)
278 : : end if
279 : 0 : call timab(48,2,tsec)
280 : : end if
281 : :
282 [ # # ]: 0 : end subroutine opernla_ylm_mv
283 : : !!***
284 : :
285 : : end module m_opernla_ylm_mv
286 : : !!***
|
