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 : !!***
|
