Line data Source code
1 : !!****m* ABINIT/m_gemm_nonlop
2 : !! NAME
3 : !! m_gemm_nonlop
4 : !!
5 : !! FUNCTION
6 : !! This module provides functions to compute the nonlocal operator by means of the BLAS GEMM
7 : !! routine. By treating ndat simultaneous wavefunctions, it is able to exploit BLAS3 routines,
8 : !! which leads to excellent CPU efficiency and OpenMP scalability.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2014-2026 ABINIT group (AL,MS)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : ! TODO list :
19 : ! Don't allocate the full nkpt structures, only those that are treated by this proc: use same init as in m_bandfft_kpt
20 : ! support more options (forces & stresses mostly)
21 : ! Support RF/other computations (only GS right now)
22 : ! handle the case where nloalg(2) < 0, ie no precomputation of ph3d
23 : ! more systematic checking of the workflow (right now, only works if init/make/gemm/destroy, no multiple makes, etc)
24 : ! Avoid allocating the complex matrix when istwfk > 1
25 : ! Merge with chebfi's invovl
26 :
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : module m_gemm_nonlop
35 :
36 : use defs_basis
37 : use m_errors
38 : use m_abicore
39 : use m_xmpi
40 : use m_xomp
41 : use m_gputk
42 : use m_abi_linalg
43 : use m_gemm_nonlop_projectors
44 :
45 : use defs_abitypes, only : MPI_type
46 : use m_opernlc_ylm, only : opernlc_ylm
47 : use m_opernlc_ylm_allwf, only : opernlc_ylm_allwf
48 : use m_opernla_gemm, only : opernla_gemm
49 : use m_opernlb_gemm, only : opernlb_gemm
50 : use m_opernld_ylm_allwf, only : opernld_ylm_allwf
51 : use m_opernld_ylm, only : opernld_ylm
52 : use m_pawcprj, only : pawcprj_type
53 : use m_geometry, only : strconv
54 : use m_kg, only : mkkpg
55 : use m_hamiltonian, only : KPRIME_H_K, K_H_KPRIME, K_H_K, KPRIME_H_KPRIME
56 :
57 : #if defined(HAVE_GPU)
58 : use m_gpu_toolbox
59 : #endif
60 :
61 : #ifdef HAVE_FC_ISO_C_BINDING
62 : use, intrinsic :: iso_c_binding, only : c_int32_t, c_int64_t, c_float, c_double, c_size_t, c_loc, c_ptr
63 : #endif
64 :
65 :
66 : implicit none
67 :
68 : private
69 :
70 : public :: gemm_nonlop
71 :
72 : ! Those routines are here to assess memory requirements
73 : public :: gemm_nonlop_ompgpu_work_mem
74 : public :: gemm_nonlop_ompgpu_static_mem
75 : !!***
76 :
77 : !----------------------------------------------------------------------
78 :
79 :
80 : contains
81 :
82 0 : function gemm_nonlop_ompgpu_work_mem(istwfk, ndat, ngrads, npw, indlmn, nattyp, ntypat, lmnmax, signs, wfoptalg) result(req_mem)
83 :
84 : integer, intent(in) :: istwfk, ndat, ngrads, npw, ntypat, lmnmax, signs, wfoptalg
85 : integer, intent(in) :: indlmn(:,:,:), nattyp(ntypat)
86 :
87 : integer :: nprojs, cplex, itypat
88 : real(dp) :: req_mem
89 :
90 : ! *************************************************************************
91 :
92 0 : cplex=2;if (istwfk>1) cplex=1
93 0 : nprojs=0
94 0 : do itypat=1,ntypat
95 0 : nprojs = nprojs + count(indlmn(3,:,itypat)>0)*nattyp(itypat)
96 : end do
97 :
98 0 : req_mem = 0
99 :
100 0 : if(cplex == 1) then
101 0 : req_mem = req_mem + dp * int(npw, c_size_t) * ndat ! temp_realvec_r
102 0 : req_mem = req_mem + dp * int(npw, c_size_t) * ndat ! temp_realvec_i
103 : end if
104 :
105 0 : req_mem = req_mem + dp * lmnmax * (lmnmax+1)/2 * ntypat ! sij_typ
106 :
107 0 : req_mem = req_mem + dp * cplex * int(nprojs, c_size_t) * int(ndat, c_size_t) ! projections
108 0 : req_mem = req_mem + dp * cplex * int(nprojs, c_size_t) * int(ndat, c_size_t) ! s_projections
109 0 : req_mem = req_mem + dp * cplex * int(nprojs, c_size_t) * int(ndat, c_size_t) ! vnl_projections
110 :
111 : ! Not in a place where vectin, vectout, svectout, enlout are allocated
112 0 : if(wfoptalg<0) then
113 0 : req_mem = req_mem + dp * cplex * int(npw, c_size_t) * int(ndat, c_size_t) ! vectin
114 0 : if(signs==2) then
115 0 : req_mem = req_mem + dp * cplex * int(npw, c_size_t) * int(ndat, c_size_t) ! vectout
116 0 : req_mem = req_mem + dp * cplex * int(npw, c_size_t) * int(ndat, c_size_t) ! svectout
117 : end if
118 0 : if(signs==1) then
119 0 : req_mem = req_mem + dp * cplex * nprojs * int(ndat, c_size_t) ! enlout (overestimate)
120 : end if
121 : end if
122 :
123 0 : if(ngrads>0) then
124 0 : req_mem = req_mem + dp * cplex * ngrads * int(nprojs, c_size_t) * int(ndat, c_size_t) ! dprojections
125 0 : if(signs==2) then
126 0 : req_mem = req_mem + dp * cplex * ngrads * int(nprojs, c_size_t) * int(ndat, c_size_t) ! s_dprojections
127 0 : req_mem = req_mem + dp * cplex * ngrads * int(nprojs, c_size_t) * int(ndat, c_size_t) ! vnl_dprojections
128 : end if
129 : end if
130 :
131 0 : end function gemm_nonlop_ompgpu_work_mem
132 :
133 : !----------------------------------------------------------------------
134 :
135 0 : function gemm_nonlop_ompgpu_static_mem(npw, indlmn, nattyp, ntypat, mpi_block_size, ngrads, use_distrib) result(req_mem)
136 :
137 : integer, intent(in) :: npw, ntypat, mpi_block_size, ngrads
138 : integer, intent(in) :: indlmn(:,:,:), nattyp(ntypat)
139 : logical, intent(in) :: use_distrib
140 :
141 : integer :: nprojs, nprojs_last_blk, itypat
142 : integer(kind=c_size_t) :: req_mem
143 :
144 : ! *************************************************************************
145 :
146 0 : nprojs = 0
147 0 : do itypat=1,ntypat
148 0 : nprojs = nprojs + count(indlmn(3,:,itypat)>0)*nattyp(itypat)
149 : end do
150 0 : nprojs_last_blk = nprojs / mpi_block_size + modulo(nprojs,mpi_block_size)
151 :
152 0 : req_mem = 0
153 :
154 0 : if(mpi_block_size>1 .and. use_distrib) then
155 0 : req_mem = req_mem + dp * 2 * int(npw, c_size_t) * int(nprojs_last_blk, c_size_t) !projs_recv
156 : #ifdef HAVE_GPU_MPI
157 : if(ngrads==0) then
158 : ! Add a suspected internal buffer for GPU-aware MPI (no derivatives)
159 : req_mem = req_mem + dp * 2 * int(npw, c_size_t) * int(nprojs_last_blk, c_size_t) !projs_recv
160 : end if
161 : #endif
162 : end if
163 : ! projs or projs_r + projs_i
164 0 : req_mem = req_mem + 2 * dp * int(npw, c_size_t) * int(nprojs_last_blk, c_size_t)
165 0 : if(ngrads>0) then
166 : ! dprojs or dprojs_r + dprojs_i
167 0 : req_mem = req_mem + 2 * dp * int(npw, c_size_t) * int(ngrads, c_size_t) * int(nprojs_last_blk, c_size_t)
168 0 : if(mpi_block_size>1 .and. use_distrib) then
169 0 : req_mem = req_mem + dp * 2 * int(npw, c_size_t) * int(ngrads, c_size_t)*int(nprojs_last_blk, c_size_t) !dprojs_recv
170 : #ifdef HAVE_GPU_MPI
171 : ! Add a suspected internal buffer for GPU-aware MPI (with derivatives)
172 : req_mem = req_mem + dp * 2 * int(npw, c_size_t) * int(ngrads, c_size_t)*int(nprojs_last_blk, c_size_t) !dprojs_recv
173 : #endif
174 : end if
175 : end if
176 :
177 0 : end function gemm_nonlop_ompgpu_static_mem
178 :
179 : !----------------------------------------------------------------------
180 :
181 : !!****f* m_gemm_nonlop/gemm_nonlop
182 : !! NAME
183 : !! gemm_nonlop
184 : !!
185 : !! FUNCTION
186 : !! Replacement of nonlop. same prototype as nonlop although not all options are implemented.
187 : !!
188 : !! INPUTS
189 : !! [gpu_option] = GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
190 : !!
191 : !! SOURCE
192 14154 : subroutine gemm_nonlop(atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,&
193 14154 : & enl,enl_ndat,enlout,ffnlin,ffnlout,gmet,gprimd,idir,indlmn,istwf_k,&
194 28308 : & kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,&
195 14154 : & mpi_enreg,natom,nattyp,ndat,ngfft,nkpgin,nkpgout,nloalg,&
196 : & nnlout,npwin,npwout,nspinor,nspinortot,ntypat,only_SO,paw_opt,&
197 28308 : & ph3din,ph3dout,signs,sij,svectout,&
198 14154 : & tim_nonlop,ucvol,useylm,vectin,vectout,&
199 14154 : & atom_proj_shift,select_k,iatom_only,typat,usepaw,&
200 14154 : & vectproj,gpu_option)
201 :
202 : !Arguments ------------------------------------
203 : !scalars
204 : integer,intent(in) :: choice,cpopt,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,idir
205 : integer,intent(in) :: istwf_k,lmnmax,matblk,mgfft,natom,ndat,nkpgin
206 : integer,intent(in) :: nkpgout,nnlout,npwin,npwout,nspinor,nspinortot,ntypat,only_SO
207 : integer,intent(in) :: paw_opt,signs,tim_nonlop,useylm,atom_proj_shift,select_k,iatom_only,usepaw
208 : integer,optional,intent(in) :: gpu_option
209 : real(dp),intent(in) :: lambda(ndat),ucvol
210 : type(MPI_type),intent(in) :: mpi_enreg
211 : !arrays
212 : integer,intent(in),target :: atindx1(natom),indlmn(6,lmnmax,ntypat),kgin(3,npwin)
213 : integer,intent(in),target :: kgout(3,npwout),nattyp(ntypat),ngfft(18),nloalg(3),typat(natom)
214 : real(dp),intent(in),ABI_CONTIGUOUS target :: enl(:,:,:,:),enl_ndat(:,:,:,:,:)
215 : real(dp),intent(in),target :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
216 : real(dp),intent(in),target :: ffnlout(npwout,dimffnlout,lmnmax,ntypat),gmet(3,3)
217 : real(dp),intent(in) :: gprimd(3,3),kptin(3),kptout(3)
218 : real(dp),intent(in),target :: kpgin(npwin,nkpgin*useylm),kpgout(npwout,nkpgout*useylm)
219 : real(dp),intent(in),target :: sij(dimenl1,ntypat*((paw_opt+1)/3))
220 : real(dp),intent(inout),target :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
221 : real(dp),intent(inout),target :: vectin(2,npwin*nspinor*ndat)
222 : real(dp),intent(inout) :: enlout(nnlout*ndat)
223 : real(dp),intent(out), target :: svectout(:,:)
224 : real(dp),intent(inout),target :: vectout(:,:)
225 : real(dp),intent(inout),optional, ABI_CONTIGUOUS target :: vectproj(:,:,:)
226 : type(pawcprj_type),intent(inout) :: cprjin(:,:)
227 :
228 : ! locals
229 : integer :: ii, idat, igrad, nprojs, ngrads, ngrads2, shift, iatom, nlmn, ierr, ibeg, iend, ikin, ikout
230 : integer :: cplex, cplex_enl, cplex_fac
231 : integer :: nnlout_test
232 : integer :: iatm, ndgxdt, ndgxdtfac, nd2gxdt, nd2gxdtfac, optder, itypat, ilmn
233 14154 : integer,allocatable :: cplex_dgxdt(:), cplex_d2gxdt(:)
234 : logical :: local_vectproj,use_enl_ndat
235 : real(dp) :: d2gxdt_dum_in(1,1,1,1,1), d2gxdt_dum_out(1,1,1,1,1),d2gxdt_dum_out2(1,1,1,1,1)
236 14154 : real(dp), allocatable :: sij_typ(:,:)
237 14154 : real(dp), ABI_CONTIGUOUS pointer :: projections(:,:,:)
238 14154 : real(dp), allocatable :: s_projections(:,:,:), vnl_projections(:,:,:)
239 14154 : real(dp), allocatable :: dprojections(:,:,:), temp_realvec_r(:), temp_realvec_i(:)
240 14154 : real(dp), allocatable, target :: s_dprojections(:,:,:), vnl_dprojections(:,:,:)
241 14154 : real(dp), allocatable, target :: d2projections(:,:,:)
242 14154 : real(dp), allocatable :: enlk(:),fnlk(:,:),ddkk(:,:),strnlk(:,:),gmet2(:,:)
243 14154 : real(dp), allocatable :: work1(:),work2(:),work3(:,:),work4(:,:),work5(:,:,:),work6(:,:,:),work7(:,:,:)
244 : integer :: idbeg,idend,idfbeg,idfend,dshift,id2beg,id2end,d2shift,dfshift,enlout_shift,ndat_enl
245 : real(dp) :: work(6)
246 : integer :: ndgxdt_stored,ishift
247 : integer :: mu0,ic,nu,mu,jc,mua,mub,nua1,nua2,nub1,nub2
248 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
249 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
250 : integer :: matblk_,natom_,ntypat_,ispden,dimenl2_,ia_beg,ia_end,dimsij,nkpgin_,nkpgout_
251 14154 : integer, ABI_CONTIGUOUS pointer :: atindx1_(:),indlmn_(:,:,:),nattyp_(:)
252 14154 : real(dp),pointer :: ffnlin_(:,:,:,:),ffnlout_(:,:,:,:)
253 14154 : real(dp),pointer :: ph3din_(:,:,:),ph3dout_(:,:,:)
254 14154 : real(dp), ABI_CONTIGUOUS pointer :: enl_(:,:,:,:,:)
255 14154 : real(dp), ABI_CONTIGUOUS pointer :: sij_(:,:)
256 14154 : real(dp), ABI_CONTIGUOUS pointer :: kpgin_(:,:),kpgout_(:,:)
257 : logical :: nld_on_gpu
258 :
259 : logical :: transfer_vectin,transfer_vectout,transfer_svectout
260 14154 : real(dp), pointer :: vectin_(:,:),vectout_(:,:),svectout_(:,:)
261 :
262 : ! *************************************************************************
263 :
264 : ! We keep the same interface as nonlop, but we don't use many of those
265 : ABI_UNUSED((/gmet/))
266 : ABI_UNUSED((/mgfft/))
267 : ABI_UNUSED((/nloalg,ngfft,only_SO,tim_nonlop/))
268 :
269 : ! Check supported options
270 14154 : if (.not.gemm_nonlop_use_gemm) then
271 0 : ABI_BUG('computation not prepared for gemm_nonlop use!')
272 : end if
273 : if ( (choice>3.and.choice/=7.and.choice/=5.and.choice/=51.and.signs==2) .or. &
274 14154 : & (choice>3.and.choice/=7.and.choice/=23.and.choice/=4.and.choice/=54.and.choice/=55.and.choice/=6.and.signs==1) .or. &
275 : & (useylm/=1) ) then
276 0 : ABI_BUG('gemm_nonlop option not supported!')
277 : end if
278 14154 : if (signs==1) then
279 2676 : nnlout_test=0
280 2676 : if (choice==1) nnlout_test=1
281 2676 : if (choice==2) nnlout_test=3*natom
282 2676 : if (choice==3) nnlout_test=6
283 2676 : if (choice==23) nnlout_test=6+3*natom
284 2676 : if (nnlout<nnlout_test) then
285 0 : ABI_BUG('wrong nnlout size!')
286 : end if
287 : end if
288 :
289 14154 : transfer_vectin=.false.; transfer_vectout=.false.; transfer_svectout=.false.
290 14154 : if(gpu_option==ABI_GPU_OPENMP) then
291 : ! Check if provided buffers are already mapped on GPU
292 : transfer_vectin=.not. xomp_target_is_present(c_loc(vectin)) &
293 : .and. ((cpopt < 2 .and. choice < 2) .or. (cpopt <= 3 .and. choice >= 2) &
294 0 : .or. (choice/=7 .and.paw_opt >=3))
295 : transfer_vectout=.not. xomp_target_is_present(c_loc(vectout)) &
296 0 : .and. (signs==2 .and. (paw_opt == 0 .or. paw_opt == 1 .or. paw_opt == 4))
297 : transfer_svectout=.not. xomp_target_is_present(c_loc(svectout)) &
298 0 : .and. (signs==2 .and. (paw_opt == 3 .or. paw_opt == 4))
299 : end if
300 :
301 14154 : ikin=1; ikout=1;
302 14154 : select case (select_k)
303 : case (K_H_K)
304 : ikin=1; ikout=1;
305 : case (K_H_KPRIME)
306 : ikin=2; ikout=1;
307 : case (KPRIME_H_K)
308 : ikin=1; ikout=2;
309 : case (KPRIME_H_KPRIME)
310 14154 : ikin=2; ikout=2;
311 : end select
312 14154 : cplex=2;if (istwf_k>1) cplex=1
313 14154 : cplex_enl=1;if (paw_opt>0) cplex_enl=2*dimenl1/(lmnmax*(lmnmax+1)) ! is enl complex?
314 14154 : cplex_fac=max(cplex,dimekbq)
315 14154 : if ((nspinortot==2.or.cplex_enl==2).and.paw_opt>0.and.choice/=7) cplex_fac=2 ! is vnl_projections complex?
316 84924 : use_enl_ndat=.false. ; if (size(enl_ndat)>0) use_enl_ndat=.true.
317 14154 : ndat_enl=1; if(use_enl_ndat) ndat_enl=ndat
318 :
319 : ! Processing one atom : set pointers to atom-specific arrays (for opernlc)
320 14154 : if(iatom_only > 0) then
321 0 : iatm=atindx1(iatom_only);itypat=typat(iatom_only)
322 0 : ia_beg=iatom_only; ia_end=iatom_only
323 0 : natom_=1 ; ntypat_=1 ; dimenl2_=1 ; matblk_=1
324 0 : ABI_MALLOC(atindx1_,(1))
325 0 : ABI_MALLOC(nattyp_,(1))
326 0 : atindx1_(1)=1 ; nattyp_(1)=1
327 0 : ABI_MALLOC(ph3din_,(2,npwin,1))
328 0 : ABI_MALLOC(ph3dout_,(2,npwout,1))
329 0 : ph3din_(:,1:npwin,1)=ph3din(:,1:npwin,iatm)
330 0 : ph3dout_(:,1:npwout,1)=ph3dout(:,1:npwout,iatm)
331 0 : ABI_MALLOC(ffnlin_,(npwin,dimffnlin,lmnmax,1))
332 0 : ABI_MALLOC(ffnlout_,(npwout,dimffnlout,lmnmax,1))
333 0 : ffnlin_(:,:,:,1)=ffnlin(:,:,:,itypat)
334 0 : ffnlout_(:,:,:,1)=ffnlout(:,:,:,itypat)
335 0 : ABI_MALLOC(indlmn_,(6,lmnmax,1))
336 0 : indlmn_(:,:,1)=indlmn(:,:,itypat)
337 0 : if (size(sij)>0) then
338 0 : dimsij=size(sij,1)
339 0 : ABI_MALLOC(sij_,(dimsij,1))
340 0 : if (size(sij,2)==ntypat) then
341 0 : sij_(:,1)=sij(:,itypat)
342 0 : else if (size(sij)>0) then
343 0 : sij_(:,1)=sij(:,1)
344 : end if
345 : end if
346 0 : if (size(enl_ndat)>0) then
347 0 : ABI_MALLOC(enl_,(size(enl_ndat,1),1,nspinor**2,ndat,size(enl_ndat,5)))
348 0 : do idat=1,ndat
349 0 : do ii=1,size(enl_ndat,5)
350 0 : do ispden=1,nspinor**2
351 0 : if (dimenl2==natom .and. usepaw==1) then
352 0 : enl_(:,1,ispden,idat,ii)=enl_ndat(:,iatom_only,ispden,idat,ii)
353 0 : else if (dimenl2==ntypat) then
354 0 : enl_(:,1,ispden,idat,ii)=enl_ndat(:,itypat,ispden,idat,ii)
355 : else
356 0 : enl_(:,1,ispden,idat,ii)=enl_ndat(:,1,ispden,idat,ii)
357 : end if
358 : end do
359 : end do
360 : end do
361 0 : else if (size(enl)>0) then
362 0 : ABI_MALLOC(enl_,(size(enl,1),1,nspinor**2,size(enl,4),1))
363 0 : do ii=1,size(enl,4)
364 0 : do ispden=1,nspinor**2
365 0 : if (dimenl2==natom .and. usepaw==1) then
366 0 : enl_(:,1,ispden,ii,1)=enl(:,iatom_only,ispden,ii)
367 0 : else if (dimenl2==ntypat) then
368 0 : enl_(:,1,ispden,ii,1)=enl(:,itypat,ispden,ii)
369 : else
370 0 : enl_(:,1,ispden,ii,1)=enl(:,1,ispden,ii)
371 : end if
372 : end do
373 : end do
374 : else
375 0 : ABI_MALLOC(enl_,(0,0,0,0,0))
376 : end if
377 :
378 : ! Usual case : all atoms are processed
379 : else
380 14154 : natom_ =natom; ntypat_=ntypat
381 14154 : ia_beg=1; ia_end=natom
382 14154 : dimenl2_=dimenl2 ; matblk_=matblk
383 14154 : atindx1_ => atindx1
384 14154 : nattyp_ => nattyp
385 14154 : ffnlin_ => ffnlin
386 14154 : ffnlout_ => ffnlout
387 14154 : enl_(1:dimenl1,1:dimenl2,1:nspinortot**2,1:dimekbq,1:1) => enl(:,:,:,:)
388 14154 : if(use_enl_ndat) then
389 0 : enl_ => enl_ndat
390 : end if
391 14154 : sij_ => sij
392 14154 : indlmn_ => indlmn
393 14154 : ph3din_ => ph3din
394 14154 : ph3dout_ => ph3dout
395 : end if
396 :
397 : ! The number of projectors used for computation may vary among
398 : ! nonlop calls, from computing on all atoms to a select one for
399 : ! some perturbations.
400 : ! In such cases, projs arrays must be recomputed
401 14154 : nprojs=0
402 37718 : do itypat=1,ntypat_
403 414430 : nprojs = nprojs + count(indlmn_(3,:,itypat)>0)*nattyp_(itypat)
404 : end do
405 :
406 14154 : if(nprojs == 0) then
407 : ! TODO check if this is correct
408 0 : if (iatom_only>0) then
409 0 : ABI_FREE(atindx1_)
410 0 : ABI_FREE(nattyp_)
411 0 : ABI_FREE(ph3din_)
412 0 : ABI_FREE(ph3dout_)
413 0 : ABI_FREE(ffnlin_)
414 0 : ABI_FREE(ffnlout_)
415 0 : ABI_FREE(enl_)
416 0 : ABI_FREE(indlmn_)
417 0 : if (size(sij) > 1) then
418 0 : ABI_FREE(sij_)
419 : end if
420 : end if
421 0 : if(signs == 1) then
422 0 : enlout=zero
423 : return
424 : end if
425 0 : if(signs == 2) then
426 0 : if(gpu_option==ABI_GPU_DISABLED) then
427 0 : vectout = zero
428 0 : if(paw_opt>0) svectout = vectin
429 : else if(gpu_option==ABI_GPU_OPENMP) then
430 : #ifdef HAVE_OPENMP_OFFLOAD
431 : if(transfer_vectout) then
432 : vectout = zero
433 : else
434 : call gpu_set_to_zero(vectout, int(2,c_size_t) * npwout * nspinor * ndat)
435 : end if
436 : if(paw_opt>0) then
437 : if(transfer_svectout .and. transfer_vectin) then
438 : svectout = vectin
439 : else if(transfer_vectin) then
440 : svectout = vectin
441 : !$OMP TARGET UPDATE TO(svectout)
442 : else if(transfer_svectout) then
443 : !$OMP TARGET UPDATE FROM(vectin)
444 : svectout = vectin
445 : else
446 : call gpu_copy(svectout, vectin, int(2,c_size_t) * npwin * nspinor * ndat)
447 : end if
448 : end if
449 : #endif
450 : end if
451 0 : return
452 : end if
453 : end if
454 :
455 : !Eventually re-compute (k+G) vectors (and related data)
456 14154 : nkpgin_=0
457 14154 : if (choice==2.or.choice==54) nkpgin_=3
458 14154 : if (signs==1) then
459 2676 : if (choice==4) nkpgin_=9
460 2676 : if (choice==3.or.choice==23.or.choice==6) nkpgin_=3
461 2676 : if (choice==55) nkpgin_=3
462 : end if
463 14154 : if (nkpgin<nkpgin_) then
464 256 : ABI_MALLOC(kpgin_,(npwin,nkpgin_))
465 64 : call mkkpg(kgin,kpgin_,kptin,nkpgin_,npwin)
466 : else
467 14090 : nkpgin_ = nkpgin
468 14090 : kpgin_ => kpgin
469 : end if
470 :
471 14154 : nkpgout_=0
472 14154 : if ((choice==2.or.choice==3.or.choice==54).and.signs==2) nkpgout_=3
473 14154 : if (nkpgout<nkpgout_) then
474 0 : ABI_MALLOC(kpgout_,(npwout,nkpgout_))
475 0 : call mkkpg(kgout,kpgout_,kptout,nkpgout_,npwout)
476 : else
477 14154 : nkpgout_ = nkpgout
478 14154 : kpgout_ => kpgout
479 : end if
480 :
481 : #ifdef HAVE_OPENMP_OFFLOAD
482 : !$OMP TARGET ENTER DATA MAP(to:kpgin_) if(nkpgin_ > 0 .and. gpu_option==ABI_GPU_OPENMP)
483 : !$OMP TARGET ENTER DATA MAP(to:kpgout_) if(nkpgout_ > 0 .and. gpu_option==ABI_GPU_OPENMP)
484 :
485 : ! Allocate and copy GPU buffers if user doesn't manage them
486 : !$OMP TARGET ENTER DATA MAP(to:vectin) IF(transfer_vectin)
487 : !$OMP TARGET ENTER DATA MAP(alloc:vectout) IF(transfer_vectout)
488 : !$OMP TARGET ENTER DATA MAP(alloc:svectout) IF(transfer_svectout)
489 :
490 : !$OMP TARGET ENTER DATA MAP(to:atindx1,indlmn) IF(gpu_option==ABI_GPU_OPENMP)
491 : if(size(enl_)>0) then
492 : !$OMP TARGET ENTER DATA MAP(to:enl_) IF(gpu_option==ABI_GPU_OPENMP)
493 : end if
494 : #endif
495 :
496 : !FIXME These seemingly useless pointers are used in BLAS operations for
497 : ! working around a bug in AOMP LLVM misreading mapped device pointers.
498 : ! I chose to generalise the workaround to avoid dupplicating each
499 : ! BLAS call specifically for handling AOMP LLVM.
500 14154 : vectin_ => vectin
501 14154 : vectout_ => vectout
502 14154 : svectout_ => svectout
503 :
504 : ! If vectproj is provided, use it for further calculations, use allocated array otherwise
505 14154 : local_vectproj=.false.
506 14154 : if(PRESENT(vectproj)) then
507 4608 : if(size(vectproj)>1) local_vectproj=.true.
508 : end if
509 1152 : if (local_vectproj) projections => vectproj
510 :
511 :
512 14154 : if(signs == 1 .and. choice > 0) then
513 18568 : enlout=zero
514 1560 : ABI_MALLOC(enlk,(ndat))
515 2824 : enlk=zero
516 2080 : ABI_MALLOC(fnlk,(3*natom,ndat))
517 48904 : fnlk=zero
518 1560 : ABI_MALLOC(ddkk,(6,ndat))
519 16648 : ddkk=zero
520 1040 : ABI_MALLOC(strnlk,(6,ndat))
521 16648 : strnlk=zero
522 : #ifdef HAVE_OPENMP_OFFLOAD
523 : !$OMP TARGET ENTER DATA MAP(to:enlout,enlk,fnlk,ddkk,strnlk) IF(gpu_option==ABI_GPU_OPENMP)
524 : #endif
525 : end if
526 :
527 14154 : if(paw_opt>=2 .and. choice > 0 .and. choice /= 7) then
528 34424 : ABI_MALLOC(sij_typ,(lmnmax*(lmnmax+1)/2,ntypat))
529 8606 : if (cplex_enl==1) then
530 22910 : do itypat=1, ntypat_
531 242696 : nlmn=count(indlmn_(3,:,itypat)>0)
532 2076314 : do ilmn=1,nlmn*(nlmn+1)/2
533 2067708 : sij_typ(ilmn,itypat)=sij_(ilmn,itypat)
534 : end do
535 : end do
536 : else
537 0 : do itypat=1, ntypat_
538 0 : nlmn=count(indlmn_(3,:,itypat)>0)
539 0 : do ilmn=1,nlmn*(nlmn+1)/2
540 0 : sij_typ(ilmn,itypat)=sij_(2*ilmn-1,itypat)
541 : end do
542 : end do
543 : end if
544 : #ifdef HAVE_OPENMP_OFFLOAD
545 : !$OMP TARGET ENTER DATA MAP(to:sij_typ) IF(gpu_option==ABI_GPU_OPENMP)
546 : #endif
547 : else
548 16644 : ABI_MALLOC(sij_typ,(1,ntypat)) ! Dummy alloc
549 : end if
550 :
551 14154 : ndgxdt = -1
552 14154 : nd2gxdt = -1
553 :
554 14154 : ndgxdtfac = 0; nd2gxdtfac = 0
555 14154 : if (choice==2) then
556 212 : if (signs==1) ndgxdt=3
557 212 : if (signs==2) ndgxdt=1
558 212 : if (signs==2) ndgxdtfac=1
559 : end if
560 14154 : if (choice==22) then
561 0 : if (signs==2) ndgxdtfac=1
562 : end if
563 14154 : if (choice==23) then
564 0 : if (signs==1) ndgxdt=9
565 : end if
566 14154 : if (choice==3) then
567 212 : if (signs==1) ndgxdt=6
568 212 : if (signs==2) ndgxdt=1
569 212 : if (signs==2) ndgxdtfac=1
570 : end if
571 14154 : if (choice==4) then
572 0 : if(signs==1) ndgxdt=3
573 0 : if(signs==1) ndgxdtfac=3
574 0 : if(signs==1) nd2gxdt=6
575 : end if
576 14154 : if (choice==5) then
577 0 : if(signs==1) ndgxdt=3
578 0 : if(signs==2) ndgxdt=1
579 0 : if(signs==2) ndgxdtfac=1
580 : end if
581 14154 : if (choice==51) then
582 0 : if(signs==1) ndgxdt=3
583 0 : if(signs==2) ndgxdt=1
584 0 : if(signs==2) ndgxdtfac=1
585 : end if
586 14154 : if (choice==54) then
587 0 : if(signs==1) ndgxdt=6
588 0 : if(signs==1) ndgxdtfac=6
589 0 : if(signs==1) nd2gxdt=9
590 0 : if(signs==2) ndgxdt=1
591 0 : if(signs==2) nd2gxdt=1
592 0 : if(signs==2) ndgxdtfac=1
593 0 : if(signs==2) nd2gxdtfac=1
594 : end if
595 14154 : if (choice==55) then
596 0 : if(signs==1) ndgxdt=9
597 0 : if(signs==1) ndgxdtfac=9
598 0 : if(signs==1) nd2gxdt=18
599 : end if
600 14154 : if (choice==6) then
601 0 : if(signs==1) ndgxdt=9
602 0 : if(signs==1) ndgxdtfac=9
603 0 : if(signs==1) nd2gxdt=54
604 : end if
605 14154 : ngrads=0; ngrads2=0
606 14154 : if(ndgxdt>0) ngrads=ndgxdt; if(ndgxdt>0) ngrads2=nd2gxdt
607 14154 : if(ndgxdt>0) then
608 424 : ABI_CHECK(ndgxdtfac<=ndgxdt,"BUG: ndgxdtfac>ndgxdt!")
609 : end if
610 14154 : optder = 0;if (ndgxdtfac>0) optder = 1
611 14154 : if (nd2gxdtfac>0) optder=2
612 42462 : ABI_MALLOC(cplex_dgxdt, (ndgxdt))
613 42462 : ABI_MALLOC(cplex_d2gxdt,(nd2gxdt))
614 16062 : cplex_dgxdt(:) = 1 ; cplex_d2gxdt(:) = 1
615 : ! When istwf_k > 1, gx derivatives can be real or pure imaginary
616 : ! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
617 : ! cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
618 14154 : if(ndgxdt > 0) then
619 424 : if (choice==5.or.choice==51) cplex_dgxdt(:) = 2
620 424 : if (choice==54.and.signs==1) cplex_dgxdt(4:6) = 2
621 : !if (choice==54.and.signs==2) cplex_dgxdt(:) = 2
622 424 : if (choice==55.and.signs==1) cplex_dgxdt(7:9) = 2
623 : end if
624 14154 : if(nd2gxdt > 0) then
625 0 : if (choice==54) cplex_d2gxdt(:) = 2
626 0 : if (choice==55.and.signs==1) cplex_d2gxdt(1:18)= 2
627 : end if
628 :
629 : ! These will store the non-local factors for vectin, svectout and vectout respectively
630 14154 : if(.not. local_vectproj) then
631 65010 : ABI_MALLOC(projections,(cplex, nprojs,nspinor*ndat))
632 : #ifdef HAVE_OPENMP_OFFLOAD
633 : !$OMP TARGET ENTER DATA MAP(alloc:projections) IF(gpu_option==ABI_GPU_OPENMP)
634 : #endif
635 : end if
636 70770 : ABI_MALLOC(s_projections,(cplex, nprojs,nspinor*ndat))
637 70770 : ABI_MALLOC(vnl_projections,(cplex_fac, nprojs,nspinor*ndat))
638 : #ifdef HAVE_OPENMP_OFFLOAD
639 : !$OMP TARGET ENTER DATA MAP(alloc:s_projections,vnl_projections) IF(gpu_option==ABI_GPU_OPENMP)
640 : #endif
641 :
642 14154 : if(gpu_option==ABI_GPU_DISABLED) then
643 8482993 : if(cpopt < 2) projections = zero
644 8694673 : s_projections = zero
645 8694673 : vnl_projections = zero
646 0 : else if(gpu_option==ABI_GPU_OPENMP) then
647 0 : if(cpopt < 2) call gpu_set_to_zero(projections, int(cplex,c_size_t)*nprojs*ndat*nspinor)
648 0 : call gpu_set_to_zero(s_projections, int(cplex,c_size_t)*nprojs*ndat*nspinor)
649 0 : call gpu_set_to_zero(vnl_projections, int(cplex_fac,c_size_t)*nprojs*ndat*nspinor)
650 : end if
651 :
652 : ! Working buffers for storing derivative
653 14154 : if (ndgxdt>0) then
654 2120 : ABI_MALLOC(dprojections,(cplex, ndgxdt*nprojs,nspinor*ndat))
655 : #ifdef HAVE_OPENMP_OFFLOAD
656 : !$OMP TARGET ENTER DATA MAP(alloc:dprojections) IF(gpu_option==ABI_GPU_OPENMP)
657 : #endif
658 424 : if(cpopt < 4) then
659 424 : if(gpu_option==ABI_GPU_DISABLED) then
660 1464640 : dprojections(:,:,:) = zero
661 0 : else if(gpu_option==ABI_GPU_OPENMP) then
662 0 : call gpu_set_to_zero(dprojections, int(cplex,c_size_t)*ndgxdt*nprojs*ndat*nspinor)
663 : end if
664 : end if
665 : else
666 41190 : ABI_MALLOC(dprojections,(1,1,ndat))
667 : end if
668 :
669 14154 : if (ndgxdtfac>0) then
670 0 : ABI_MALLOC(s_dprojections,(cplex, ndgxdtfac*nprojs,nspinor*ndat))
671 0 : ABI_MALLOC(vnl_dprojections,(cplex_fac, ndgxdtfac*nprojs,nspinor*ndat))
672 : #ifdef HAVE_OPENMP_OFFLOAD
673 : !$OMP TARGET ENTER DATA MAP(alloc:s_dprojections,vnl_dprojections) IF(gpu_option==ABI_GPU_OPENMP)
674 : #endif
675 0 : if(gpu_option==ABI_GPU_DISABLED) then
676 0 : s_dprojections(:,:,:) = zero
677 0 : vnl_dprojections(:,:,:) = zero
678 0 : else if(gpu_option==ABI_GPU_OPENMP) then
679 0 : call gpu_set_to_zero(s_dprojections, int(cplex,c_size_t)*ndgxdtfac*nprojs*ndat*nspinor)
680 0 : call gpu_set_to_zero(vnl_dprojections, int(cplex_fac,c_size_t)*ndgxdtfac*nprojs*ndat*nspinor)
681 : end if
682 : else
683 42462 : ABI_MALLOC(s_dprojections,(1,1,ndat))
684 28308 : ABI_MALLOC(vnl_dprojections,(1,1,ndat))
685 : end if
686 :
687 : ! Working buffers for storing 2nd-derivative
688 14154 : if (nd2gxdt>0) then
689 0 : ABI_MALLOC(d2projections,(cplex, nd2gxdt*nprojs, nspinor*ndat))
690 : #ifdef HAVE_OPENMP_OFFLOAD
691 : !$OMP TARGET ENTER DATA MAP(alloc:d2projections) IF(gpu_option==ABI_GPU_OPENMP)
692 : #endif
693 0 : if(cpopt < 4) then
694 0 : if(gpu_option==ABI_GPU_DISABLED) then
695 0 : d2projections(:,:,:) = zero
696 0 : else if(gpu_option==ABI_GPU_OPENMP) then
697 0 : call gpu_set_to_zero(d2projections, int(cplex,c_size_t)*nd2gxdt*nprojs*ndat*nspinor)
698 : end if
699 : end if
700 : else
701 42462 : ABI_MALLOC(d2projections,(1, 1, ndat))
702 : end if
703 :
704 :
705 : ! determine precisely when temp_realvec_r~i needs to be allocated
706 : ! to factorize allocate (resp. deallocate) at the begining (resp. at the end) of subroutine
707 : ! to avoid multiple allocate/deallocate that can be costly
708 14154 : if (cplex /= 2) then
709 : if ( (cpopt < 2) .or. &
710 4111 : & (paw_opt == 3 .or. paw_opt == 4) .or. &
711 : & (paw_opt == 0 .or. paw_opt == 1 .or. paw_opt == 4)) then
712 12333 : ABI_MALLOC(temp_realvec_r,(MAX(npwout,npwin)*nspinor*ndat))
713 8222 : ABI_MALLOC(temp_realvec_i,(MAX(npwout,npwin)*nspinor*ndat))
714 : #ifdef HAVE_OPENMP_OFFLOAD
715 : !$OMP TARGET ENTER DATA MAP(alloc:temp_realvec_r,temp_realvec_i) IF(gpu_option==ABI_GPU_OPENMP)
716 : #endif
717 : end if
718 : end if
719 :
720 14154 : if(cpopt >= 2) then
721 : ! retrieve from cprjin
722 864 : if(.not. local_vectproj .and. cpopt/=3) then
723 : !$OMP PARALLEL DO PRIVATE(shift,idat,iatom,nlmn)
724 2016 : do idat=1, ndat*nspinor
725 : shift = 0
726 5472 : do iatom = ia_beg, ia_end
727 3456 : nlmn = cprjin(iatom, idat)%nlmn
728 86400 : projections(1:cplex, shift+1:shift+nlmn, idat) = cprjin(iatom, idat)%cp(1:cplex, 1:nlmn)
729 5184 : shift = shift + nlmn
730 : end do
731 : end do
732 : #ifdef HAVE_OPENMP_OFFLOAD
733 : !$OMP TARGET UPDATE TO(projections) IF(gpu_option==ABI_GPU_OPENMP)
734 : #endif
735 : end if
736 864 : if(cpopt==4.and.allocated(dprojections)) then
737 0 : ABI_CHECK(cprjin(1,1)%ncpgr>=ndgxdt,"cprjin%ncpgr not correct! (1)")
738 0 : ndgxdt_stored = cprjin(1,1)%ncpgr
739 0 : ishift=0
740 0 : if (((choice==2).or.(choice==3)).and.(ndgxdt_stored>ndgxdt).and.(signs==2)) ishift=idir-ndgxdt
741 0 : if ((choice==2).and.(ndgxdt_stored==9).and.(signs==2)) ishift=ishift+6
742 0 : if (choice==2.and.(ndgxdt_stored>ndgxdt).and.(signs==1)) ishift=ndgxdt_stored-ndgxdt
743 : !$OMP PARALLEL DO PRIVATE(shift,idat,iatom,igrad,nlmn)
744 0 : do idat=1, ndat*nspinor
745 : shift = 0
746 0 : do iatom = ia_beg, ia_end
747 0 : nlmn = cprjin(iatom, idat)%nlmn
748 0 : do ilmn=1,nlmn
749 0 : do igrad=1,ndgxdt
750 : dprojections(1:cplex, shift + igrad, idat) = &
751 0 : cprjin(iatom, idat)%dcp(1:cplex,igrad+ishift,ilmn)
752 : end do
753 0 : shift = shift + ndgxdt
754 : end do
755 : end do
756 : end do
757 : #ifdef HAVE_OPENMP_OFFLOAD
758 : !$OMP TARGET UPDATE TO(dprojections) IF(gpu_option==ABI_GPU_OPENMP)
759 : #endif
760 : end if
761 : end if ! cpopt
762 :
763 14154 : if(cpopt<=1.or.(cpopt<=3.and.(choice==2.or.choice==3.or.choice==5.or.choice==51.or.choice==23.or.choice==54.or.choice==55.or.choice==4))) then
764 :
765 : call opernla_gemm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnlin,&
766 : & d2projections,dprojections,ffnlin,projections,&
767 : & idir,indlmn,istwf_k,kpgin_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nkpgin_,&
768 : & npwin,nspinor,ph3din,signs,ucvol,ndat,ntypat,lmnmax,nattyp,(ikin==2),&
769 : & iatom_only,atom_proj_shift,cpopt,&
770 : & nprojs,&
771 : & vectin,&
772 : & temp_realvec_r,temp_realvec_i,&
773 13290 : & gpu_option,gemm_nonlop_is_distributed)
774 :
775 13290 : if(cpopt >= 0) then
776 : ! store in cprjin
777 2780 : if(.not. local_vectproj .and. cpopt/=3) then
778 : #ifdef HAVE_OPENMP_OFFLOAD
779 : !$OMP TARGET UPDATE FROM(projections) IF(gpu_option==ABI_GPU_OPENMP)
780 : #endif
781 : !$OMP PARALLEL DO PRIVATE(shift,idat,iatom,nlmn)
782 5492 : do idat=1, ndat*nspinor
783 : shift = 0
784 12068 : do iatom = ia_beg, ia_end
785 6576 : nlmn = cprjin(iatom, idat)%nlmn
786 242656 : cprjin(iatom, idat)%cp(1:cplex, 1:nlmn) = projections(1:cplex, shift+1:shift+nlmn, idat)
787 9864 : shift = shift + nlmn
788 : end do
789 : end do
790 : end if
791 2780 : if(cpopt==1 .or. cpopt==3) then
792 0 : ABI_CHECK(cprjin(1,1)%ncpgr>=ndgxdt,"cprjin%ncpgr not correct! (2)")
793 : #ifdef HAVE_OPENMP_OFFLOAD
794 : !$OMP TARGET UPDATE FROM(dprojections) IF(gpu_option==ABI_GPU_OPENMP)
795 : #endif
796 : !$OMP PARALLEL DO PRIVATE(shift,idat,iatom,igrad,nlmn)
797 0 : do idat=1, ndat*nspinor
798 : shift = 0
799 0 : do iatom = ia_beg, ia_end
800 0 : nlmn = cprjin(iatom, idat)%nlmn
801 0 : do ilmn=1,nlmn
802 0 : do igrad=1,ndgxdt
803 : cprjin(iatom, idat)%dcp(1:cplex,igrad,ilmn) = &
804 0 : & dprojections(1:cplex, shift + igrad, idat)
805 : end do
806 0 : shift = shift + ndgxdt
807 : end do
808 : end do
809 : end do
810 : end if
811 : end if ! cpopt >= 0
812 : end if ! cpopt >= 2
813 :
814 14154 : if(choice > 0) then
815 :
816 11422 : if(choice /= 7) then
817 : ! opernlc
818 10846 : iatm = 0
819 10846 : shift = 0; dshift = 0; dfshift = 0; d2shift = 0
820 29342 : do itypat=1, ntypat_
821 319464 : nlmn=count(indlmn_(3,:,itypat)>0)
822 :
823 18496 : ibeg = shift+1
824 18496 : iend = shift+nattyp_(itypat)*nlmn
825 :
826 18496 : idbeg = dshift+1
827 18496 : idend = dshift+nattyp_(itypat)*nlmn*ngrads
828 :
829 18496 : idfbeg = dshift+1
830 18496 : idfend = dshift+nattyp_(itypat)*nlmn*ndgxdtfac
831 :
832 18496 : id2beg = d2shift+1
833 18496 : id2end = d2shift+nattyp_(itypat)*nlmn*ngrads2
834 :
835 : call opernlc_ylm_allwf(atindx1_,cplex,cplex_dgxdt,cplex_d2gxdt,&
836 : & cplex_enl,cplex_fac,&
837 : & dprojections,&
838 : & vnl_dprojections,&
839 : & s_dprojections,&
840 : & d2projections,d2gxdt_dum_out,d2gxdt_dum_out2,&
841 : & dimenl1,dimenl2_,dimekbq,enl_,&
842 : & projections,&
843 : & vnl_projections,&
844 : & s_projections,&
845 : & iatm,indlmn_(:,:,itypat),itypat,lambda,mpi_enreg,natom_,&
846 : & ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,&
847 : & nattyp_(itypat),nlmn,nspinor,nspinortot,optder,paw_opt,sij_typ(:,itypat),&
848 18496 : & ndat,ibeg-1,iend,nprojs,ndat_enl,gpu_option)
849 :
850 18496 : shift = shift + nattyp_(itypat)*nlmn
851 18496 : dshift = dshift + nattyp_(itypat)*nlmn*ngrads
852 18496 : dfshift = dshift + nattyp_(itypat)*nlmn*ndgxdtfac
853 18496 : d2shift = d2shift + nattyp_(itypat)*nlmn*ngrads2
854 29342 : iatm = iatm+nattyp_(itypat)
855 : end do
856 : else
857 576 : if(gpu_option==ABI_GPU_DISABLED) then
858 128160 : s_projections = projections
859 0 : else if(gpu_option==ABI_GPU_OPENMP) then
860 0 : call gpu_copy(s_projections, projections, int(cplex,c_size_t) * nprojs * nspinor * ndat)
861 : end if
862 : end if ! choice /= 7
863 :
864 : ! opernlb
865 11422 : if(signs==2) then
866 :
867 : call opernlb_gemm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
868 : & d2gxdt_dum_in,d2gxdt_dum_out,&
869 : & vnl_dprojections,s_dprojections,dimffnlout,ffnlout,&
870 : & vnl_projections,s_projections,&
871 : & idir,indlmn,kpgout_,matblk,istwf_k,&
872 : & nd2gxdt,nd2gxdtfac,ndgxdt,ndgxdtfac,&
873 : & nkpgout_,npwout,nspinor,signs,ucvol,ndat,&
874 : & ntypat,lmnmax,nattyp,(ikout==2),iatom_only,atom_proj_shift,&
875 : & paw_opt,ph3dout,&
876 : & nprojs,&
877 : & vectin_,vectout_,svectout_,&
878 : & temp_realvec_r,temp_realvec_i,&
879 10902 : & gpu_option,gemm_nonlop_is_distributed)
880 : end if
881 :
882 : ! opernld
883 11422 : if(signs==1) then
884 520 : nld_on_gpu = .false.
885 520 : if(choice==1 .or. choice==2 .or. choice==3 .or. choice==23 .or. choice==4 .or. choice==54 .or. choice==55 .or. choice==6) then
886 520 : if(gpu_option==ABI_GPU_OPENMP) nld_on_gpu = .true.
887 : call opernld_ylm_allwf(choice,cplex,cplex_fac,ddkk,&
888 : & dprojections,vnl_dprojections,s_dprojections,d2projections,&
889 : & enlk,enlout,fnlk,projections,vnl_projections,s_projections,&
890 : & natom,ndat,nd2gxdt,ndgxdt,&
891 : & ndgxdtfac,indlmn_,ntypat_,lmnmax,nprojs,nnlout,nspinor,paw_opt,&
892 520 : & strnlk,nattyp_,gpu_option)
893 : else
894 0 : shift=0; dshift=0; dfshift = 0; d2shift = 0; iatm=1
895 : #ifdef HAVE_OPENMP_OFFLOAD
896 : !$OMP TARGET UPDATE FROM(dprojections,vnl_dprojections,s_dprojections) IF(gpu_option==ABI_GPU_OPENMP)
897 : !$OMP TARGET UPDATE FROM(d2projections) IF(gpu_option==ABI_GPU_OPENMP)
898 : !$OMP TARGET UPDATE FROM(projections,vnl_projections,s_projections) IF(gpu_option==ABI_GPU_OPENMP)
899 : #endif
900 0 : do itypat=1, ntypat_
901 0 : nlmn=count(indlmn_(3,:,itypat)>0)
902 :
903 0 : ibeg = shift+1
904 0 : iend = shift+nattyp_(itypat)*nlmn
905 :
906 0 : idbeg = dshift+1
907 0 : idend = dshift+nattyp_(itypat)*nlmn*ngrads
908 :
909 0 : idfbeg = dshift+1
910 0 : idfend = dshift+nattyp_(itypat)*nlmn*ndgxdtfac
911 :
912 0 : id2beg = d2shift+1
913 0 : id2end = d2shift+nattyp_(itypat)*nlmn*ngrads2
914 :
915 0 : do idat=1,ndat
916 : call opernld_ylm (choice,cplex,cplex_fac,ddkk(:,idat),&
917 : & dprojections (:, idbeg:idend, 1+nspinor*(idat-1):nspinor*idat),&
918 : & vnl_dprojections(:, idfbeg:idfend, 1+nspinor*(idat-1):nspinor*idat),&
919 : & s_dprojections (:, idfbeg:idfend, 1+nspinor*(idat-1):nspinor*idat),&
920 : & d2projections (:, id2beg:id2end, 1+nspinor*(idat-1):nspinor*idat),&
921 : & enlk(idat),enlout(nnlout*(idat-1)+1:nnlout*idat),fnlk(:,idat),&
922 : & projections (:, ibeg:iend, 1+nspinor*(idat-1):nspinor*idat),&
923 : & vnl_projections(:, ibeg:iend, 1+nspinor*(idat-1):nspinor*idat),&
924 : & s_projections (:, ibeg:iend, 1+nspinor*(idat-1):nspinor*idat),&
925 : & iatm,natom_,1,nd2gxdt,ndgxdt,ndgxdtfac,&
926 0 : & nattyp_(itypat),nlmn,nnlout,nspinor,paw_opt,strnlk(:,idat))
927 : end do
928 :
929 0 : shift = shift + nattyp_(itypat)*nlmn
930 0 : dshift = dshift + nattyp_(itypat)*nlmn*ngrads
931 0 : dfshift = dshift + nattyp_(itypat)*nlmn*ndgxdtfac
932 0 : d2shift = d2shift + nattyp_(itypat)*nlmn*ngrads2
933 0 : iatm = iatm+nattyp_(itypat)
934 : end do
935 : end if
936 :
937 : #ifdef HAVE_OPENMP_OFFLOAD
938 : !$OMP TARGET UPDATE FROM(enlout) if(nld_on_gpu)
939 : #endif
940 :
941 : ! Reduction in case of parallelism
942 520 : if (mpi_enreg%paral_spinor==1) then
943 0 : if (size(enlout)>0) then
944 0 : call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr)
945 : end if
946 0 : if (choice==3.or.choice==23) then
947 : #ifdef HAVE_OPENMP_OFFLOAD
948 : !$OMP TARGET UPDATE FROM(enlk) if(nld_on_gpu)
949 : #endif
950 0 : call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr)
951 : end if
952 0 : if (choice==55) then
953 0 : call xmpi_sum(ddkk,mpi_enreg%comm_spinor,ierr)
954 : end if
955 : end if
956 :
957 : !Need sometimes gmet
958 520 : if ((signs==1.and.paw_opt<=3).and. &
959 : & (choice==5 .or.choice==51.or.choice==52.or.choice==53.or.&
960 : & choice==54.or.choice==55)) then
961 0 : ABI_MALLOC(gmet2,(3,3))
962 520 : gmet2 = MATMUL(TRANSPOSE(gprimd),gprimd)
963 : end if
964 :
965 : !Coordinate transformations
966 :
967 : ! Derivatives wrt strain
968 : ! - Convert from reduced to cartesian coordinates
969 : ! - Substract volume contribution
970 520 : if ((choice==3.or.choice==23).and.paw_opt<=3) then
971 : #ifdef HAVE_OPENMP_OFFLOAD
972 : !$OMP TARGET UPDATE FROM(enlk) if(nld_on_gpu)
973 : #endif
974 980 : do idat=1,ndat
975 768 : enlout_shift=(idat-1)*nnlout
976 768 : call strconv(enlout(enlout_shift+1:enlout_shift+6),gprimd,work)
977 3072 : enlout(enlout_shift+1:enlout_shift+3)=(work(1:3)-enlk(idat))
978 3284 : enlout(enlout_shift+4:enlout_shift+6)= work(4:6)
979 : end do
980 : end if
981 :
982 : !2nd derivative wrt to k wave vector and atomic position (effective charges):
983 : ! - convert from cartesian to reduced coordinates
984 520 : if (choice==54.and.signs==1.and.paw_opt<=3) then
985 0 : ABI_MALLOC(work1,(3))
986 0 : ABI_MALLOC(work2,(3))
987 0 : do idat=1,ndat
988 0 : mu0=0 ! Shift to be applied in enlout array
989 0 : enlout_shift=(idat-1)*nnlout
990 0 : do mu=1,3*natom
991 : ! First, real part
992 0 : work1(1)=enlout(enlout_shift+mu0+1);work1(2)=enlout(enlout_shift+mu0+3);work1(3)=enlout(enlout_shift+mu0+5)
993 0 : work2(:)=gmet2(:,1)*work1(1)+gmet2(:,2)*work1(2)+gmet2(:,3)*work1(3)
994 0 : enlout(enlout_shift+mu0+1)=work2(1);enlout(enlout_shift+mu0+3)=work2(2);enlout(enlout_shift+mu0+5)=work2(3)
995 : ! Then imaginary part
996 0 : work1(1)=enlout(enlout_shift+mu0+2);work1(2)=enlout(enlout_shift+mu0+4);work1(3)=enlout(enlout_shift+mu0+6)
997 0 : work2(:)=gmet2(:,1)*work1(1)+gmet2(:,2)*work1(2)+gmet2(:,3)*work1(3)
998 0 : enlout(enlout_shift+mu0+2)=work2(1);enlout(enlout_shift+mu0+4)=work2(2);enlout(enlout_shift+mu0+6)=work2(3)
999 0 : mu0=mu0+6
1000 : end do
1001 : end do !idat
1002 0 : ABI_FREE(work1)
1003 0 : ABI_FREE(work2)
1004 : end if
1005 :
1006 : !2nd derivative wrt to k wave vector and strain (piezoelectric tensor):
1007 : ! - convert from cartesian to reduced coordinates (k point)
1008 : ! - convert from reduced to cartesian coordinates (strain)
1009 : ! - substract volume contribution
1010 : ! - symetrize strain components
1011 520 : if (choice==55.and.signs==1.and.paw_opt<=3) then
1012 : #ifdef HAVE_OPENMP_OFFLOAD
1013 : !$OMP TARGET UPDATE FROM(ddkk) if(nld_on_gpu)
1014 : #endif
1015 0 : ABI_MALLOC(work3,(2,3))
1016 0 : ABI_MALLOC(work4,(2,3))
1017 0 : ABI_MALLOC(work5,(2,3,6))
1018 0 : ABI_MALLOC(work7,(2,3,6))
1019 0 : ABI_MALLOC(work6,(2,3,3))
1020 0 : do idat=1,ndat
1021 0 : enlout_shift=(idat-1)*nnlout
1022 0 : do ic=1,3 ! gamma
1023 0 : work5=zero
1024 0 : do jc=1,3 ! nu
1025 0 : do ii=1,3 ! lambda
1026 0 : mu=(gamma(jc,ii)-1)*3+1
1027 : work5(1,jc,ii)=gmet2(ic,1)*enlout(enlout_shift+2*mu-1)+gmet2(ic,2)*enlout(enlout_shift+2*mu+1) &
1028 0 : & +gmet2(ic,3)*enlout(enlout_shift+2*mu+3)
1029 : work5(2,jc,ii)=gmet2(ic,1)*enlout(enlout_shift+2*mu )+gmet2(ic,2)*enlout(enlout_shift+2*mu+2) &
1030 0 : & +gmet2(ic,3)*enlout(enlout_shift+2*mu+4)
1031 : end do
1032 : end do
1033 0 : work6=zero
1034 0 : do jc=1,3 ! nu
1035 0 : do ii=1,3 ! beta
1036 : work6(1:cplex,ii,jc)=gprimd(ii,1)*work5(1:cplex,jc,1)+gprimd(ii,2)*work5(1:cplex,jc,2) &
1037 0 : & +gprimd(ii,3)*work5(1:cplex,jc,3)
1038 : end do
1039 : end do
1040 0 : do jc=1,3 ! alpha
1041 0 : do ii=1,3 ! beta
1042 0 : mu=gamma(jc,ii)
1043 : work7(1:cplex,ic,mu)=gprimd(jc,1)*work6(1:cplex,ii,1)+gprimd(jc,2)*work6(1:cplex,ii,2) &
1044 0 : & +gprimd(jc,3)*work6(1:cplex,ii,3)
1045 : end do
1046 : end do
1047 : end do ! gamma
1048 :
1049 0 : do ii=1,3 ! alpha
1050 : work3(1,ii)=gprimd(ii,1)*ddkk(2*1-1,idat)+gprimd(ii,2)*ddkk(2*2-1,idat) &
1051 0 : & +gprimd(ii,3)*ddkk(2*3-1,idat)
1052 : work3(2,ii)=gprimd(ii,1)*ddkk(2*1 ,idat)+gprimd(ii,2)*ddkk(2*2 ,idat) &
1053 0 : & +gprimd(ii,3)*ddkk(2*3 ,idat)
1054 : end do
1055 0 : do ii=1,3 ! gamma
1056 : work4(1,ii)=gmet2(ii,1)*ddkk(2*1-1,idat)+gmet2(ii,2)*ddkk(2*2-1,idat) &
1057 0 : & +gmet2(ii,3)*ddkk(2*3-1,idat)
1058 : work4(2,ii)=gmet2(ii,1)*ddkk(2*1 ,idat)+gmet2(ii,2)*ddkk(2*2 ,idat) &
1059 0 : & +gmet2(ii,3)*ddkk(2*3 ,idat)
1060 : end do
1061 :
1062 0 : do mu=1,6
1063 0 : ii=alpha(mu) ! alpha
1064 0 : ic=beta(mu) ! beta
1065 0 : do jc=1,3 ! gamma
1066 : work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-half &
1067 0 : & *(gprimd(ic,jc)*work3(1:cplex,ii)+gprimd(ii,jc)*work3(1:cplex,ic))
1068 0 : if (ii==ic) work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-work4(1:cplex,jc)
1069 : end do
1070 : end do
1071 0 : do mu=1,6 ! alpha,beta
1072 0 : do nu=1,3 ! gamma
1073 0 : mu0=3*(mu-1)+nu
1074 0 : enlout(enlout_shift+2*mu0-1)=work7(1,nu,mu)
1075 0 : enlout(enlout_shift+2*mu0 )=work7(2,nu,mu)
1076 : end do
1077 : end do
1078 : end do !idat
1079 0 : ABI_FREE(work3)
1080 0 : ABI_FREE(work4)
1081 0 : ABI_FREE(work5)
1082 0 : ABI_FREE(work6)
1083 0 : ABI_FREE(work7)
1084 : end if
1085 :
1086 :
1087 : !2nd derivative wrt to 2 strains (elastic tensor):
1088 : ! - convert from reduced to cartesian coordinates
1089 : ! - substract volume contribution
1090 520 : if (choice==6.and.signs==1.and.paw_opt<=3) then
1091 : #ifdef HAVE_OPENMP_OFFLOAD
1092 : !$OMP TARGET UPDATE FROM(enlk,strnlk,fnlk) if(nld_on_gpu)
1093 : #endif
1094 0 : ABI_MALLOC(work1,(6))
1095 0 : ABI_MALLOC(work2,(6))
1096 0 : ABI_MALLOC(work3,(6+3*natom,6))
1097 0 : do idat=1,ndat
1098 0 : mu0=(idat-1)*nnlout ! Shift to be applied in enlout array
1099 0 : work3(:,:)=reshape(enlout(mu0+1:mu0+6*(6+3*natom)),(/6+3*natom,6/))
1100 0 : do mu=1,6
1101 0 : call strconv(work3(1:6,mu),gprimd,work3(1:6,mu))
1102 : end do
1103 0 : do mu=1,6+3*natom
1104 0 : work1(1:6)=work3(mu,1:6)
1105 0 : call strconv(work1,gprimd,work2)
1106 0 : work3(mu,1:6)=work2(1:6)
1107 : end do
1108 0 : enlout(mu0+1:mu0+6*(6+3*natom))=reshape(work3(:,:),(/6*(6+3*natom)/))
1109 0 : call strconv(strnlk(:,idat),gprimd,strnlk(:,idat))
1110 0 : do mub=1,6
1111 0 : nub1=alpha(mub);nub2=beta(mub)
1112 0 : do mua=1,6
1113 0 : mu=mu0+mua+(3*natom+6)*(mub-1)
1114 0 : nua1=alpha(mua);nua2=beta(mua)
1115 0 : if (mua<=3.and.mub<=3) enlout(mu)=enlout(mu)+enlk(idat)
1116 0 : if (mua<=3) enlout(mu)=enlout(mu)-strnlk(mub,idat)
1117 0 : if (mub<=3) enlout(mu)=enlout(mu)-strnlk(mua,idat)
1118 0 : if (nub1==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub2),idat)
1119 0 : if (nub2==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub1),idat)
1120 0 : if (nub1==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub2),idat)
1121 0 : if (nub2==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub1),idat)
1122 : end do
1123 0 : if (mub<=3) then
1124 0 : do nua1=1,natom
1125 0 : nua2=3*(nua1-1);mu=mu0+nua2+6+(3*natom+6)*(mub-1)
1126 0 : enlout(mu+1:mu+3)=enlout(mu+1:mu+3)-fnlk(nua2+1:nua2+3,idat)
1127 : end do
1128 : end if
1129 : end do
1130 : end do
1131 0 : ABI_FREE(work1)
1132 0 : ABI_FREE(work2)
1133 0 : ABI_FREE(work3)
1134 : end if
1135 :
1136 : end if !opernld
1137 :
1138 : end if ! choice>0
1139 :
1140 : #ifdef HAVE_OPENMP_OFFLOAD
1141 : if(gpu_option==ABI_GPU_OPENMP) then
1142 : ! Retrieve and release allocated buffers
1143 : !$OMP TARGET EXIT DATA MAP(delete:vectin) IF(transfer_vectin)
1144 : !$OMP TARGET EXIT DATA MAP(from:vectout) IF(transfer_vectout)
1145 : !$OMP TARGET EXIT DATA MAP(from:svectout) IF(transfer_svectout)
1146 :
1147 : !$OMP TARGET EXIT DATA MAP(delete:s_projections,vnl_projections)
1148 : !$OMP TARGET EXIT DATA MAP(delete:projections) IF(.not. local_vectproj)
1149 : !$OMP TARGET EXIT DATA MAP(delete:dprojections) IF(ndgxdt>0)
1150 : !$OMP TARGET EXIT DATA MAP(delete:s_dprojections) IF(ndgxdtfac>0)
1151 : !$OMP TARGET EXIT DATA MAP(delete:vnl_dprojections) IF(ndgxdtfac>0)
1152 : !$OMP TARGET EXIT DATA MAP(delete:d2projections) IF(nd2gxdt>0)
1153 :
1154 : if (cplex /= 2) then
1155 : if ( (cpopt < 2) .or. &
1156 : & (paw_opt == 3 .or. paw_opt == 4) .or. &
1157 : & (paw_opt == 0 .or. paw_opt == 1 .or. paw_opt == 4)) then
1158 : !$OMP TARGET EXIT DATA MAP(delete:temp_realvec_r,temp_realvec_i)
1159 : end if
1160 : end if
1161 :
1162 : if(size(enl_)>0) then
1163 : !$OMP TARGET EXIT DATA MAP(delete:enl_)
1164 : end if
1165 :
1166 : if(paw_opt>=2 .and. choice > 0 .and. choice /= 7) then
1167 : !$OMP TARGET EXIT DATA MAP(delete:sij_typ) IF(gpu_option==ABI_GPU_OPENMP)
1168 : end if
1169 :
1170 : !$OMP TARGET EXIT DATA MAP(delete:kpgin_) IF(nkpgin_ > 0)
1171 : !$OMP TARGET EXIT DATA MAP(delete:kpgout_) IF(nkpgout_ > 0)
1172 :
1173 : !$OMP TARGET EXIT DATA MAP(delete:enlk,fnlk,strnlk,ddkk,enlout) IF(signs == 1 .and. choice > 0)
1174 : !$OMP TARGET EXIT DATA MAP(delete:atindx1,indlmn)
1175 : end if
1176 : #endif
1177 :
1178 : ! Release memory
1179 :
1180 14154 : if (iatom_only>0) then
1181 0 : ABI_FREE(atindx1_)
1182 0 : ABI_FREE(nattyp_)
1183 0 : ABI_FREE(ph3din_)
1184 0 : ABI_FREE(ph3dout_)
1185 0 : ABI_FREE(ffnlin_)
1186 0 : ABI_FREE(ffnlout_)
1187 0 : ABI_FREE(enl_)
1188 0 : ABI_FREE(indlmn_)
1189 0 : if (size(sij) > 1) then
1190 0 : ABI_FREE(sij_)
1191 : end if
1192 : end if
1193 :
1194 14154 : if(signs == 1 .and. choice > 0) then
1195 : #ifdef HAVE_OPENMP_OFFLOAD
1196 : #endif
1197 520 : ABI_FREE(enlk)
1198 520 : ABI_FREE(fnlk)
1199 520 : ABI_FREE(strnlk)
1200 520 : ABI_FREE(ddkk)
1201 : end if
1202 :
1203 14154 : if (nkpgin<nkpgin_) then
1204 64 : ABI_FREE(kpgin_)
1205 : end if
1206 14154 : if (nkpgout<nkpgout_) then
1207 0 : ABI_FREE(kpgout_)
1208 : end if
1209 :
1210 14154 : if (allocated(gmet2)) then
1211 0 : ABI_FREE(gmet2)
1212 : end if
1213 :
1214 14154 : if (allocated(sij_typ)) then
1215 14154 : ABI_FREE(sij_typ)
1216 : end if
1217 :
1218 14154 : ABI_FREE(cplex_dgxdt)
1219 14154 : ABI_FREE(cplex_d2gxdt)
1220 :
1221 14154 : if(.not. local_vectproj) then
1222 13002 : ABI_FREE(projections)
1223 : end if
1224 14154 : ABI_FREE(s_projections)
1225 14154 : ABI_FREE(vnl_projections)
1226 14154 : if (allocated(dprojections)) then
1227 14154 : ABI_FREE(dprojections)
1228 : end if
1229 14154 : if (allocated(s_dprojections)) then
1230 14154 : ABI_FREE(s_dprojections)
1231 : end if
1232 14154 : if (allocated(vnl_dprojections)) then
1233 14154 : ABI_FREE(vnl_dprojections)
1234 : end if
1235 14154 : if (allocated(d2projections)) then
1236 14154 : ABI_FREE(d2projections)
1237 : end if
1238 14154 : if (allocated(temp_realvec_r)) then
1239 4111 : ABI_FREE(temp_realvec_r)
1240 4111 : ABI_FREE(temp_realvec_i)
1241 : end if
1242 :
1243 28308 : end subroutine gemm_nonlop
1244 : !***
1245 :
1246 : end module m_gemm_nonlop
1247 : !!***
|
