Branch data Line data Source code
1 : : !!****m* ABINIT/m_oper
2 : : !! NAME
3 : : !! m_oper
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
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 : : !! INPUTS
14 : : !!
15 : : !! OUTPUT
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : ! nvtx related macro definition
27 : : #include "nvtx_macros.h"
28 : :
29 : : MODULE m_oper
30 : :
31 : : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
32 : : use defs_basis
33 : : use m_abicore
34 : : use m_errors
35 : : use m_xomp
36 : : use m_gputk
37 : : use m_abi_linalg
38 : :
39 : : !use m_abi_linalg, only : abi_xgemm
40 : : use m_hide_lapack, only : xginv
41 : : use m_matlu, only : copy_matlu,destroy_matlu,diff_matlu,identity_matlu,init_matlu, &
42 : : & inverse_matlu,matlu_type,print_matlu,prod_matlu,trace_matlu,zero_matlu
43 : : use m_paw_dmft, only : mpi_distrib_dmft_type,paw_dmft_type
44 : : use m_xmpi, only : xmpi_allgatherv,xmpi_gatherv,xmpi_sum,xmpi_sum_master
45 : :
46 : : #ifdef HAVE_GPU_MARKERS
47 : : use m_nvtx_data
48 : : #endif
49 : :
50 : : implicit none
51 : :
52 : : private
53 : :
54 : : public :: init_oper
55 : : public :: init_oper_ndat
56 : : public :: diff_oper
57 : : public :: destroy_oper
58 : : public :: print_oper
59 : : public :: inverse_oper
60 : : public :: downfold_oper
61 : : public :: identity_oper
62 : : public :: copy_oper
63 : : public :: copy_oper_from_ndat
64 : : public :: copy_oper_to_ndat
65 : : public :: trace_oper
66 : : public :: upfold_oper
67 : : public :: prod_oper
68 : : public :: trace_prod_oper
69 : : public :: gather_oper
70 : : public :: gather_oper_ks
71 : : !!***
72 : :
73 : : !!****t* m_oper/oper_type
74 : : !! NAME
75 : : !! oper_type
76 : : !!
77 : : !! FUNCTION
78 : : !! This structured datatype contains the necessary data
79 : : !!
80 : : !! SOURCE
81 : :
82 : : type, public :: oper_type ! for each atom
83 : :
84 : : ! integer :: maxlpawu ! Number of correlated atoms
85 : : !
86 : : ! integer :: mband
87 : : ! ! Number of bands
88 : :
89 : : ! Wether ks and matlu are stored on GPU
90 : : integer :: gpu_option
91 : :
92 : : integer :: ndat
93 : :
94 : : integer :: has_operks
95 : : ! Is the operator allocated in the KS basis ?
96 : :
97 : : integer :: has_opermatlu
98 : : ! Is the operator allocated in the local basis ?
99 : : !
100 : : integer :: mbandc
101 : : ! Total number of correlated bands
102 : :
103 : : integer :: natom
104 : : ! Number of atoms
105 : :
106 : : integer :: nkpt
107 : : ! Number of k-point in the IBZ.
108 : : !
109 : : integer :: nspinor
110 : : ! Number of spinors
111 : : !
112 : : integer :: nsppol
113 : : ! Number of spin polarizations
114 : :
115 : : integer :: paral
116 : : ! =1 if the operator has been memory-parallelized over kpt, 0 otherwise
117 : :
118 : : integer :: shiftk
119 : : ! Shift to get the physical kpt index (when the operator is memory-parallelized over kpt)
120 : :
121 : : !character(len=12) :: whichoper
122 : : ! describe the type of operator computed (DFT, DMFT, KS..)
123 : :
124 : : ! ! Polarisation
125 : : type(matlu_type), allocatable :: matlu(:)
126 : : ! Local projection on correlated orbitals
127 : :
128 : : complex(dp), allocatable :: ks(:,:,:,:)
129 : : ! In the KS basis (mbandc,mbandc,nkpt,nsppol)
130 : :
131 : : real(dp), ABI_CONTIGUOUS pointer :: wtk(:) => null()
132 : : ! Weights for each kpt
133 : :
134 : : end type oper_type
135 : : !!***
136 : :
137 : : !----------------------------------------------------------------------
138 : :
139 : :
140 : : CONTAINS !========================================================================================
141 : : !!***
142 : :
143 : : !!****f* m_oper/init_oper
144 : : !! NAME
145 : : !! init_oper
146 : : !!
147 : : !! FUNCTION
148 : : !! Allocate variables used in type oper_type.
149 : : !!
150 : : !! INPUTS
151 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
152 : : !! nkpt = number of k-pts
153 : : !! wtk = weights for each k-pt
154 : : !! shiftk = shift for the kpt index
155 : : !! opt_ksloc = 1: initialize in KS space only
156 : : !! = 2: initialize in local space only
157 : : !! = 3: initialize in both KS and local space
158 : : !!
159 : : !! OUTPUTS
160 : : !! oper <type(oper_type)>= operator
161 : : !!
162 : : !! SOURCE
163 : :
164 : 4 : subroutine init_oper(paw_dmft,oper,nkpt,wtk,shiftk,opt_ksloc)
165 : :
166 : : !Arguments ------------------------------------
167 : : integer, optional, intent(in) :: nkpt,opt_ksloc,shiftk
168 : : type(paw_dmft_type), intent(in) :: paw_dmft
169 : : type(oper_type), intent(inout) :: oper
170 : : real(dp), target, optional :: wtk(paw_dmft%nkpt)
171 : : !Local variables ------------------------------------
172 : : integer :: optksloc
173 : : !************************************************************************
174 : :
175 : : DBG_ENTER("COLL")
176 : :
177 : 231029 : optksloc = 3
178 [ + + ]: 231029 : if (present(opt_ksloc)) optksloc = opt_ksloc
179 : :
180 : : !if(optksloc/=3) then
181 : : ! FIXME: empty line!
182 : : !endif
183 : :
184 : 231029 : oper%has_operks = 0
185 : 231029 : oper%has_opermatlu = 0
186 : 231029 : oper%gpu_option = ABI_GPU_DISABLED
187 : :
188 : : ! ===================
189 : : ! Integers
190 : : ! ===================
191 : 231029 : oper%mbandc = paw_dmft%mbandc
192 : 231029 : oper%natom = paw_dmft%natom
193 : 231029 : oper%nspinor = paw_dmft%nspinor
194 : 231029 : oper%nsppol = paw_dmft%nsppol
195 : 231029 : oper%paral = 0
196 : 231029 : oper%shiftk = 0
197 : 231029 : oper%ndat = 1
198 : :
199 [ - + ]: 231029 : if (present(shiftk)) oper%shiftk = shiftk
200 : :
201 : 231029 : oper%nkpt = paw_dmft%nkpt
202 [ + + ]: 231029 : if (present(nkpt)) oper%nkpt = nkpt
203 : :
204 [ + - + + ]: 231029 : if (present(shiftk) .or. oper%nkpt /= paw_dmft%nkpt) oper%paral = 1
205 : :
206 : : ! allocate(oper%wtk(oper%nkpt))
207 [ + + ]: 231029 : if (present(wtk)) then
208 : 4 : oper%wtk => wtk(:)
209 : : else
210 : 231025 : oper%wtk => paw_dmft%wtk(:)
211 : : end if ! present(wtk)
212 : :
213 : : ! ===================
214 : : ! KS variables
215 : : ! ===================
216 [ + + ]: 231029 : if (optksloc == 1 .or. optksloc == 3) then
217 : :
218 [ - + + - : 75978 : ABI_MALLOC(oper%ks,(oper%mbandc,oper%mbandc,oper%nkpt,oper%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
219 : 12663 : oper%has_operks = 1
220 [ + + + + : 13107306 : oper%ks(:,:,:,:) = czero
+ + + + ]
221 : :
222 : : end if ! optksloc=1 or optksloc=3
223 : :
224 : : ! ===================
225 : : ! matlu variables
226 : : ! ===================
227 [ + - ]: 231029 : if (optksloc == 2 .or. optksloc == 3) then
228 [ - + + - : 1109170 : ABI_MALLOC(oper%matlu,(oper%natom))
+ - + - -
+ - + - +
+ + + - ]
229 : 231029 : oper%has_opermatlu = 1
230 : 231029 : call init_matlu(oper%natom,oper%nspinor,oper%nsppol,paw_dmft%lpawu(:),oper%matlu(:))
231 : : end if ! optksloc=2 or optksloc=3
232 : :
233 : : DBG_EXIT("COLL")
234 : :
235 [ + + ]: 231029 : end subroutine init_oper
236 : : !!***
237 : :
238 : : !!****f* m_oper/init_oper_ndat
239 : : !! NAME
240 : : !! init_oper_ndat
241 : : !!
242 : : !! FUNCTION
243 : : !! Allocate variables used in type oper_type.
244 : : !!
245 : : !! INPUTS
246 : : !!
247 : : !! OUTPUTS
248 : : !! oper = operator of type oper_type
249 : : !!
250 : : !! SOURCE
251 : :
252 : 0 : subroutine init_oper_ndat(paw_dmft,oper,ndat,nkpt,wtk,shiftk,opt_ksloc,gpu_option)
253 : :
254 : : use m_matlu, only : init_matlu
255 : : use m_paw_dmft, only : paw_dmft_type
256 : :
257 : : !Arguments ------------------------------------
258 : : integer, optional, intent(in) :: nkpt,opt_ksloc,shiftk,gpu_option
259 : : integer, intent(in) :: ndat
260 : : type(paw_dmft_type), intent(in) :: paw_dmft
261 : : type(oper_type), target, intent(inout) :: oper
262 : : real(dp), target, optional :: wtk(paw_dmft%nkpt)
263 : : !Local variables ------------------------------------
264 : : integer :: optksloc,ndat_,l_gpu_option
265 : : !************************************************************************
266 : :
267 : : DBG_ENTER("COLL")
268 : :
269 : 0 : optksloc = 3
270 [ # # ]: 0 : if (present(opt_ksloc)) optksloc = opt_ksloc
271 [ # # ]: 0 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option=gpu_option
272 : :
273 : : !if(optksloc/=3) then
274 : : ! FIXME: empty line!
275 : : !endif
276 : 0 : ndat_=ndat;
277 : :
278 : 0 : oper%gpu_option = l_gpu_option
279 : 0 : oper%has_operks = 0
280 : 0 : oper%has_opermatlu = 0
281 : :
282 : : ! ===================
283 : : ! Integers
284 : : ! ===================
285 : 0 : oper%mbandc = paw_dmft%mbandc
286 : 0 : oper%natom = paw_dmft%natom
287 : 0 : oper%nspinor = paw_dmft%nspinor
288 : 0 : oper%nsppol = paw_dmft%nsppol
289 : 0 : oper%paral = 0
290 : 0 : oper%shiftk = 0
291 : 0 : oper%ndat = ndat_
292 : :
293 [ # # ]: 0 : if (present(shiftk)) oper%shiftk = shiftk
294 : :
295 : 0 : oper%nkpt = paw_dmft%nkpt
296 [ # # ]: 0 : if (present(nkpt)) oper%nkpt = nkpt
297 : :
298 [ # # # # ]: 0 : if (present(shiftk) .or. oper%nkpt /= paw_dmft%nkpt) oper%paral = 1
299 : :
300 : : ! allocate(oper%wtk(oper%nkpt))
301 [ # # ]: 0 : if (present(wtk)) then
302 : 0 : oper%wtk => wtk(:)
303 : : else
304 : 0 : oper%wtk => paw_dmft%wtk(:)
305 : : end if ! present(wtk)
306 : :
307 : : ! ===================
308 : : ! KS variables
309 : : ! ===================
310 [ # # ]: 0 : if (optksloc == 1 .or. optksloc == 3) then
311 : :
312 [ # # # # : 0 : ABI_MALLOC(oper%ks,(oper%mbandc,oper%mbandc*ndat_,oper%nkpt,oper%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
313 : 0 : oper%has_operks = 1
314 : : #ifdef HAVE_OPENMP_OFFLOAD
315 : : if(l_gpu_option==ABI_GPU_OPENMP) call oper_ompgpu_map_alloc(oper%ks)
316 : : #endif
317 [ # # ]: 0 : if(gpu_option==ABI_GPU_OPENMP) then
318 : 0 : call gpu_set_to_zero_complex(oper%ks, int(oper%nsppol,c_size_t)*ndat_*oper%mbandc*oper%mbandc*oper%nkpt)
319 : : else
320 [ # # # # : 0 : oper%ks(:,:,:,:) = czero
# # # # ]
321 : : end if
322 : :
323 : : end if ! optksloc=1 or optksloc=3
324 : :
325 : : ! ===================
326 : : ! matlu variables
327 : : ! ===================
328 [ # # ]: 0 : if (optksloc == 2 .or. optksloc == 3) then
329 [ # # # # : 0 : ABI_MALLOC(oper%matlu,(oper%natom))
# # # # #
# # # # #
# # ]
330 : 0 : oper%has_opermatlu = 1
331 : 0 : call init_matlu(oper%natom,oper%nspinor,oper%nsppol*ndat_,paw_dmft%lpawu(:),oper%matlu(:),gpu_option=l_gpu_option)
332 : : end if ! optksloc=2 or optksloc=3
333 : :
334 : : DBG_EXIT("COLL")
335 : :
336 : : #if defined HAVE_OPENMP_OFFLOAD
337 : : contains
338 : : ! ks is located in an array of structs,
339 : : ! which syntax is loosely supported in OpenMP directives
340 : : subroutine oper_ompgpu_map_alloc(ks)
341 : : complex(dp),intent(inout) :: ks(:,:,:,:)
342 : : !$OMP TARGET ENTER DATA MAP(alloc:ks)
343 : : end subroutine oper_ompgpu_map_alloc
344 : : #endif
345 [ # # ]: 0 : end subroutine init_oper_ndat
346 : : !!***
347 : :
348 : : !!****f* m_oper/destroy_oper
349 : : !! NAME
350 : : !! destroy_oper
351 : : !!
352 : : !! FUNCTION
353 : : !! deallocate oper
354 : : !!
355 : : !! INPUTS
356 : : !! oper <type(oper_type)>= operator
357 : : !!
358 : : !! OUTPUT
359 : : !!
360 : : !! SOURCE
361 : :
362 : 231029 : subroutine destroy_oper(oper)
363 : :
364 : : !Arguments ------------------------------------
365 : : type(oper_type), target, intent(inout) :: oper
366 : : !Local variables-------------------------------
367 : : !! *********************************************************************
368 : :
369 : : DBG_ENTER("COLL")
370 : :
371 [ + - ]: 231029 : if (oper%has_opermatlu == 1) then
372 : 231029 : call destroy_matlu(oper%matlu(:),oper%natom)
373 : : !else
374 : : ! message = " Operator is not defined to be used in destroy_oper"
375 : : ! ABI_ERROR(message)
376 : : end if ! has_opermatlu=1
377 : :
378 [ + - ]: 231029 : if (allocated(oper%matlu)) then
379 [ + + - + : 647112 : ABI_FREE(oper%matlu)
- + ]
380 : 231029 : oper%has_opermatlu = 0
381 : : end if
382 : :
383 [ + + ]: 231029 : if (allocated(oper%ks)) then
384 : : #ifdef HAVE_OPENMP_OFFLOAD
385 : : if(oper%gpu_option==ABI_GPU_OPENMP) call oper_ompgpu_map_delete(oper%ks)
386 : : #endif
387 : 12663 : ABI_FREE(oper%ks)
388 : 12663 : oper%has_operks = 0
389 : : end if
390 : :
391 : 231029 : oper%wtk => null()
392 : : ! no deallocation for wtk: wtk is an explicit pointer
393 : :
394 : : DBG_EXIT("COLL")
395 : :
396 : : #if defined HAVE_OPENMP_OFFLOAD
397 : : contains
398 : : ! ks is located in an array of structs,
399 : : ! which syntax is loosely supported in OpenMP directives
400 : : subroutine oper_ompgpu_map_delete(ks)
401 : : complex(dp),intent(inout) :: ks(:,:,:,:)
402 : : !$OMP TARGET EXIT DATA MAP(delete:ks)
403 : : end subroutine oper_ompgpu_map_delete
404 : : #endif
405 : :
406 : 0 : end subroutine destroy_oper
407 : : !!***
408 : :
409 : : !!****f* m_oper/copy_oper
410 : : !! NAME
411 : : !! copy_oper
412 : : !!
413 : : !! FUNCTION
414 : : !! Copy oper1 into oper2
415 : : !!
416 : : !! INPUTS
417 : : !! oper1 <type(oper_type)>= operator
418 : : !!
419 : : !! OUTPUT
420 : : !! oper2 <type(oper_type)>= operator
421 : : !!
422 : : !! SOURCE
423 : :
424 : 18063 : subroutine copy_oper(oper1,oper2)
425 : :
426 : : !Arguments ------------------------------------
427 : : type(oper_type), intent(in) :: oper1
428 : : type(oper_type), intent(inout) :: oper2 !vz_i
429 : : !Local variables-------------------------------
430 : : ! *********************************************************************
431 : :
432 : : DBG_ENTER("COLL")
433 : :
434 [ + - + - ]: 18063 : if (oper1%has_opermatlu == 1 .and. oper2%has_opermatlu == 1) then
435 : 18063 : call copy_matlu(oper1%matlu(:),oper2%matlu(:),oper1%natom)
436 : : end if
437 : :
438 [ + + + - ]: 18063 : if (oper1%has_operks == 1 .and. oper2%has_operks == 1) &
439 [ + + + + : 55773 : & oper2%ks(:,:,:,:) = oper1%ks(:,:,:,:)
+ + + + ]
440 : :
441 : : DBG_EXIT("COLL")
442 : :
443 : 18063 : end subroutine copy_oper
444 : : !!***
445 : :
446 : : !!****f* m_oper/copy_oper_from_ndat
447 : : !! NAME
448 : : !! copy_oper_from_ndat
449 : : !!
450 : : !! FUNCTION
451 : : !!
452 : : !! INPUTS
453 : : !!
454 : : !! OUTPUT
455 : : !!
456 : : !! SOURCE
457 : :
458 : 0 : subroutine copy_oper_from_ndat(oper1,oper2,ndat,nw,proct,me_freq,copy_ks)
459 : :
460 : : use defs_basis
461 : : use m_matlu, only : copy_matlu_from_ndat
462 : : use m_errors
463 : :
464 : : !Arguments ------------------------------------
465 : : !type
466 : : integer,intent(in) :: nw,ndat,me_freq
467 : : logical,intent(in) :: copy_ks
468 : : integer,intent(in) :: proct(nw)
469 : : type(oper_type),target,intent(in) :: oper1
470 : : type(oper_type),intent(inout) :: oper2(nw) !vz_i
471 : :
472 : : !oper variables-------------------------------
473 : : integer :: ikpt, isppol, idat, iw, iatom, mbandc
474 : 0 : complex(dp), ABI_CONTIGUOUS pointer :: mat(:,:,:)
475 : : ! *********************************************************************
476 : : DBG_ENTER("COLL")
477 : 0 : ABI_CHECK(oper1%ndat==ndat, "Bad value for ndat!")
478 : 0 : mbandc=oper1%mbandc
479 : : if(oper1%has_opermatlu==1 .and. oper1%gpu_option==ABI_GPU_OPENMP) then
480 : : do iatom=1,oper1%natom
481 : : if (oper1%matlu(iatom)%lpawu == -1) cycle
482 : : mat => oper1%matlu(iatom)%mat ! array of structs in OpenMP loosely supported
483 : : #ifdef HAVE_OPENMP_OFFLOAD
484 : : !$OMP TARGET UPDATE FROM(mat)
485 : : #endif
486 : : end do
487 : : end if
488 : 0 : idat=1
489 [ # # ]: 0 : do iw=1,nw
490 [ # # ]: 0 : if (proct(iw) /= me_freq) cycle
491 [ # # # # ]: 0 : if(oper1%has_opermatlu==1.and.oper2(iw)%has_opermatlu==1) then
492 : 0 : call copy_matlu_from_ndat(oper1%matlu,oper2(iw)%matlu,oper1%natom,ndat,idat)
493 : 0 : idat=idat+1
494 : : endif
495 : : enddo
496 : :
497 [ # # # # ]: 0 : if(allocated(oper1%ks) .and. copy_ks) then
498 : : #ifdef HAVE_OPENMP_OFFLOAD
499 : : !$OMP TARGET UPDATE FROM(oper1%ks) IF(oper1%gpu_option==ABI_GPU_OPENMP)
500 : : #endif
501 : 0 : idat=1
502 [ # # ]: 0 : do iw=1,nw
503 [ # # ]: 0 : if (proct(iw) /= me_freq) cycle
504 [ # # ]: 0 : do isppol=1,oper1%nsppol
505 [ # # ]: 0 : do ikpt=1,oper1%nkpt
506 [ # # # # ]: 0 : oper2(iw)%ks(:,:,ikpt,isppol)=oper1%ks(:,1+(idat-1)*mbandc:idat*mbandc,ikpt,isppol)
507 : : enddo
508 : : enddo
509 : 0 : idat=idat+1
510 : : enddo
511 : : endif
512 : :
513 : : DBG_EXIT("COLL")
514 [ # # ]: 0 : end subroutine copy_oper_from_ndat
515 : : !!***
516 : :
517 : : !!****f* m_oper/copy_oper_to_ndat
518 : : !! NAME
519 : : !! copy_oper_to_ndat
520 : : !!
521 : : !! FUNCTION
522 : : !!
523 : : !! INPUTS
524 : : !!
525 : : !! OUTPUT
526 : : !!
527 : : !! SOURCE
528 : :
529 : 0 : subroutine copy_oper_to_ndat(oper1,oper2,ndat,nw,proct,me_freq,copy_ks)
530 : :
531 : 0 : use defs_basis
532 : : use m_matlu, only : copy_matlu_to_ndat
533 : : use m_errors
534 : :
535 : : !Arguments ------------------------------------
536 : : !type
537 : : integer,intent(in) :: nw,ndat,me_freq
538 : : logical,intent(in) :: copy_ks
539 : : integer,intent(in) :: proct(nw)
540 : : type(oper_type),intent(in) :: oper1(nw)
541 : : type(oper_type),target,intent(inout) :: oper2 !vz_i
542 : :
543 : : !oper variables-------------------------------
544 : : integer :: ikpt, isppol, idat, iw, iatom, mbandc
545 : 0 : complex(dp), ABI_CONTIGUOUS pointer :: mat(:,:,:)
546 : : ! *********************************************************************
547 : : DBG_ENTER("COLL")
548 : 0 : ABI_CHECK(oper2%ndat==ndat, "Bad value for ndat!")
549 [ # # ]: 0 : ABI_CHECK(oper2%mbandc==oper1(1)%mbandc, "Bad value for mbandc!")
550 : 0 : mbandc=oper2%mbandc
551 : 0 : idat=1
552 [ # # ]: 0 : do iw=1,nw
553 [ # # ]: 0 : if (proct(iw) /= me_freq) cycle
554 [ # # # # ]: 0 : if(oper1(iw)%has_opermatlu==1.and.oper2%has_opermatlu==1) then
555 : 0 : call copy_matlu_to_ndat(oper1(iw)%matlu,oper2%matlu,oper2%natom,ndat,idat)
556 : 0 : idat=idat+1
557 : : endif
558 : : enddo
559 : : if(oper2%has_opermatlu==1 .and. oper2%gpu_option==ABI_GPU_OPENMP) then
560 : : do iatom=1,oper2%natom
561 : : if (oper2%matlu(iatom)%lpawu == -1) cycle
562 : : mat => oper2%matlu(iatom)%mat ! array of structs in OpenMP loosely supported
563 : : #ifdef HAVE_OPENMP_OFFLOAD
564 : : !$OMP TARGET UPDATE TO(mat)
565 : : #endif
566 : : end do
567 : : end if
568 : :
569 [ # # # # ]: 0 : if(allocated(oper2%ks) .and. copy_ks) then
570 [ # # ]: 0 : ABI_CHECK(size(oper2%ks,dim=2) == mbandc*ndat, "well?")
571 [ # # ]: 0 : ABI_CHECK(size(oper2%ks,dim=1) == mbandc, "uh?")
572 : 0 : idat=1
573 [ # # ]: 0 : do iw=1,nw
574 [ # # ]: 0 : if (proct(iw) /= me_freq) cycle
575 [ # # ]: 0 : do isppol=1,oper2%nsppol
576 [ # # ]: 0 : do ikpt=1,oper2%nkpt
577 [ # # # # ]: 0 : oper2%ks(:,1+(idat-1)*mbandc:idat*mbandc,ikpt,isppol)=oper1(iw)%ks(:,:,ikpt,isppol)
578 : : enddo
579 : : enddo
580 : 0 : idat=idat+1
581 : : enddo
582 : : #ifdef HAVE_OPENMP_OFFLOAD
583 : : !$OMP TARGET UPDATE TO(oper2%ks) IF(oper2%gpu_option==ABI_GPU_OPENMP)
584 : : #endif
585 : : endif
586 : :
587 : : DBG_EXIT("COLL")
588 [ # # ]: 0 : end subroutine copy_oper_to_ndat
589 : : !!***
590 : :
591 : : !!****f* m_oper/print_oper
592 : : !! NAME
593 : : !! print_oper
594 : : !!
595 : : !! FUNCTION
596 : : !!
597 : : !! INPUTS
598 : : !! oper <type(oper_type)>= operator
599 : : !! option= < 5: write diagonal part of KS occupation matrix
600 : : !! > 8: write all elements of KS occup. matrix.
601 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
602 : : !! prtopt= in local space: print option for print_matlu
603 : : !! in KS space: only prints if abs(prtopt)>=3
604 : : !! print off-diagonal elements if abs(prtopt)>=4
605 : : !!
606 : : !! OUTPUT
607 : : !!
608 : : !! SOURCE
609 : :
610 : 1261 : subroutine print_oper(oper,option,paw_dmft,prtopt)
611 : :
612 : : !Arguments ------------------------------------
613 : : type(paw_dmft_type), intent(in) :: paw_dmft
614 : : type(oper_type), intent(in) :: oper
615 : : integer, intent(in) :: option,prtopt
616 : : !Local variables-------------------------------
617 : : integer :: ib,ib1,iband1,iband2,ikpt,isppol,mbandc,nkpt,nkptr
618 : : character(len=50000) :: message
619 : : logical :: ximag
620 : : real(dp) :: maximag
621 : : ! *********************************************************************
622 : :
623 : : DBG_ENTER("COLL")
624 : :
625 [ + - ]: 1261 : if (oper%has_opermatlu == 1) then
626 : 1261 : write(message,'(2a)') ch10,' = In the atomic basis'
627 : 1261 : call wrtout(std_out,message,'COLL')
628 : 1261 : call print_matlu(oper%matlu(:),oper%natom,prtopt)
629 : : end if ! has_opermatlu=1
630 : :
631 [ + + ]: 1261 : if (oper%has_operks == 1) then
632 : 582 : write(message,'(2a)') ch10,' = In the Kohn-Sham basis'
633 : 582 : call wrtout(std_out,message,'COLL')
634 : :
635 : : !todo_ba complete print_out
636 : 582 : mbandc = oper%mbandc
637 : 582 : iband1 = 1
638 : 582 : iband2 = mbandc
639 : 582 : maximag = zero
640 : 582 : nkpt = oper%nkpt
641 : 582 : ximag = .false.
642 : : ! do ib=1,oper%mbandc
643 : : ! if(-(paw_dmft%eigen_dft(1,1,ib)+paw_dmft%fermie).ge.0.3) iband1=ib
644 : : ! if( (paw_dmft%eigen_dft(1,1,ib)-paw_dmft%fermie).le.0.3) iband2=ib
645 : : ! enddo
646 : :
647 [ + - + + ]: 582 : if (abs(prtopt) >= 3 .and. ((option < 5) .or. (option > 8))) then
648 : : ! write(message,'(x,a,a,i4,2x,a)') ch10,' -KS states'
649 : : ! call wrtout(std_out,message,'COLL')
650 : 186 : nkptr = min(nkpt,4)
651 [ + + ]: 448 : do isppol=1,paw_dmft%nsppol
652 : 262 : write(message,'(a,3x,a,1x,i1)') ch10,"--isppol--",isppol
653 : 262 : call wrtout(std_out,message,'COLL')
654 : 262 : write(message,'(2a)') ch10,&
655 : 524 : & " - (in the following only the values for the correlated bands and the first k-points are printed)"
656 : 262 : call wrtout(std_out,message,'COLL')
657 [ + + ]: 1356 : do ikpt=1,nkptr
658 : 908 : write(message,'(2a,i4,2x,f14.5,a)') ch10,&
659 : 1816 : & " -k-pt--",ikpt,oper%wtk(ikpt),"(<-weight(k-pt))"
660 : 908 : call wrtout(std_out,message,'COLL')
661 [ + + ]: 908 : if (option < 5) then
662 : 454 : write(message,'(19x,a,6x,a)') "Eigenvalues","Occupations"
663 : 454 : call wrtout(std_out,message,'COLL')
664 : : else if (abs(prtopt) >= 4 .or. option > 8) then
665 [ + - + + ]: 5878 : write(message,'(a,10x,2000(i5,12x))') ch10,(paw_dmft%include_bands(ib),ib=iband1,iband2)
666 : 454 : call wrtout(std_out,message,'COLL')
667 : : end if ! option
668 [ + + ]: 12018 : do ib=1,mbandc
669 [ + + ]: 10848 : if (option < 5) then
670 [ - + ]: 5424 : if (abs(aimag(oper%ks(ib,ib,ikpt,isppol))) >= tol10) then
671 : 0 : write(message,'(a,i5,e14.5,3x,e14.5,3x,e21.14)') " -iband--",paw_dmft%include_bands(ib),&
672 : 0 : & paw_dmft%eigen_dft(ib,ikpt,isppol),oper%ks(ib,ib,ikpt,isppol)
673 : : else
674 : 5424 : write(message,'(a,i5,e14.5,3x,e14.5)') " -iband--",paw_dmft%include_bands(ib),&
675 : 10848 : & paw_dmft%eigen_dft(ib,ikpt,isppol),dble(oper%ks(ib,ib,ikpt,isppol))
676 : : end if ! imaginary part
677 : 5424 : call wrtout(std_out,message,'COLL')
678 : : end if ! option<5
679 [ + - + + ]: 10848 : if (abs(prtopt) >= 4 .or. option > 8 .and. ib >= iband1 .and. ib <= iband2) then
680 : :
681 : 86784 : write(message,'(i5,1x,2000(2f7.3,3x))') paw_dmft%include_bands(ib),(dble(oper%ks(ib,ib1,ikpt,isppol)), &
682 [ + - + + ]: 92208 : & aimag(oper%ks(ib,ib1,ikpt,isppol)),ib1=iband1,iband2)
683 : 5424 : call wrtout(std_out,message,'COLL')
684 : :
685 : : ! to write imaginary part
686 : : ! write(message, '(1000(2f9.3,2x))') &
687 : : !& (real(oper%ks(isppol,ikpt,ib,ib1)),imag(oper%ks(isppol,ikpt,ib,ib1)),ib1=iband1,iband2)
688 : : ! call wrtout(std_out,message,'COLL')
689 : : end if ! prtopt>=20
690 [ + - ]: 11756 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then ! no sense to perform this check on off-diagonal elements
691 [ # # ]: 0 : if (abs(aimag(oper%ks(ib,ib,ikpt,isppol))) > max(tol10,maximag)) then
692 : 0 : ximag = .true.
693 : 0 : maximag = aimag(oper%ks(ib,ib,ikpt,isppol))
694 : : end if
695 : : else
696 [ + + ]: 184416 : do ib1=1,mbandc
697 [ + + ]: 184416 : if (abs(aimag(oper%ks(ib1,ib,ikpt,isppol))) > max(tol10,maximag)) then
698 : 544 : ximag = .true.
699 : 544 : maximag = aimag(oper%ks(ib1,ib,ikpt,isppol))
700 : : end if
701 : : end do ! ib1
702 : : end if
703 : : end do ! ib
704 : : end do ! ikpt
705 : : end do ! isppol
706 : : else
707 : 396 : write(message,'(5x,a,i10,a)') '(not written)'
708 : 396 : call wrtout(std_out,message,'COLL')
709 : : end if ! abs(prtopt)>=3 and (option<5 or option>8)
710 [ + + ]: 186 : if (ximag) then
711 : 72 : write(message,'(3a,e12.4,a)') "Occupations are imaginary !",ch10,&
712 : 144 : & " Maximal value is ",maximag,ch10
713 : 72 : ABI_WARNING(message)
714 : : end if ! ximag
715 [ + + - + ]: 679 : else if (abs(prtopt) >= 3 .and. ((option < 5) .or. (option > 8))) then
716 : 0 : write(message, '(2a)') ch10," Prb with options and has_operks in print_oper"
717 : 0 : call wrtout(std_out,message,'COLL')
718 : : end if ! if oper%has_operks
719 : : ! write(message, '(2a)') ch10," end print_oper"
720 : : ! call wrtout(std_out,message,'COLL')
721 : :
722 : : DBG_EXIT("COLL")
723 : :
724 : 0 : end subroutine print_oper
725 : : !!***
726 : :
727 : : !!****f* m_oper/inverse_oper
728 : : !! NAME
729 : : !! inverse_oper
730 : : !!
731 : : !! FUNCTION
732 : : !! Compute the inverse of the operator either in the KS space or in the
733 : : !! correlated subspace.
734 : : !!
735 : : !! INPUTS
736 : : !! oper <type(oper_type)>= operator
737 : : !! option=1 do inversion in KS band space
738 : : !! =2 do inversion in local space
739 : : !! =3 do both
740 : : !! procb(ikpt)=for kpt parallelization; gives the rank (in the kpt communicator) of the CPU handling each ikpt
741 : : !! iproc=rank of the current process in the kpt communicator
742 : : !!
743 : : !! OUTPUT
744 : : !! oper <type(oper_type)>= operator inverted
745 : : !!
746 : : !! SOURCE
747 : :
748 : 982154 : subroutine inverse_oper(oper,option,procb,iproc,gpu_option)
749 : :
750 : : !Arguments ------------------------------------
751 : : integer, intent(in) :: option
752 : : type(oper_type), target, intent(inout) :: oper
753 : : integer, optional, intent(in) :: iproc,gpu_option
754 : : integer, optional, intent(in) :: procb(oper%nkpt)
755 : : !Local variables-------------------------------
756 : : integer :: ikpt,isppol,idat,paral,mbandc
757 : : integer :: l_gpu_option
758 : : !integer :: blk,iatom,ib
759 : 982154 : complex(dp), ABI_CONTIGUOUS pointer :: ks(:,:,:,:)
760 : : #ifdef HAVE_OPENMP_OFFLOAD
761 : : complex(dp), allocatable :: work(:,:)
762 : : complex(dp), ABI_CONTIGUOUS pointer :: mat(:,:,:)
763 : : #endif
764 : : !todo_ba: prb with gwpc here: necessary for matcginv but should be dp
765 : : ! *********************************************************************
766 : :
767 : : DBG_ENTER("COLL")
768 : : ABI_NVTX_START_RANGE(NVTX_DMFT_INVERSE_OPER)
769 : :
770 : 0 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option=gpu_option
771 : 982154 : paral = 0
772 : 982154 : ks => oper%ks
773 : 982154 : mbandc = oper%mbandc
774 [ + + + + ]: 982154 : if (present(procb) .and. present(iproc) .and. oper%paral == 0) paral = 1
775 : :
776 : : !if (((option == 1 .or. option == 3) .and. (oper%has_operks == 0)) .or. &
777 : : ! & ((option == 2 .or. option == 3) .and. (oper%has_opermatlu == 0))) then
778 : : ! message = " Options are not coherent with definitions of this operator"
779 : : ! ABI_ERROR(message)
780 : : !end if
781 : :
782 [ + + ]: 982154 : if (option == 2 .or. option == 3) then
783 : 32249 : call inverse_matlu(oper%matlu(:),oper%natom)
784 : : end if
785 : :
786 : : #ifdef HAVE_OPENMP_OFFLOAD
787 : : !$OMP TARGET ENTER DATA MAP(alloc:ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
788 : : !$OMP TARGET UPDATE TO(ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
789 : : #endif
790 [ + + ]: 982154 : if (option == 1 .or. option == 3) then
791 [ + + ]: 2546255 : do isppol=1,oper%nsppol
792 [ + + ]: 11011555 : do ikpt=1,oper%nkpt
793 [ + + ]: 8465300 : if (paral == 1) then
794 [ + + ]: 132000 : if (procb(ikpt) /= iproc) cycle
795 : : end if
796 [ + - ]: 9962650 : if(l_gpu_option==ABI_GPU_DISABLED) then
797 [ + + ]: 16732600 : do idat=1,oper%ndat
798 : : ! write(std_out,*) "isppol,ikpt",isppol,ikpt,m
799 : : ! write(std_out,*) "isppol,ikpt",matrix
800 : : !call matcginv_dpc(matrix,oper%mbandc,oper%mbandc)
801 : 16732600 : call xginv(oper%ks(:,1+(idat-1)*oper%mbandc:idat*oper%mbandc,ikpt,isppol),oper%mbandc)
802 : : end do ! idat
803 : : else if(l_gpu_option==ABI_GPU_OPENMP) then
804 : : #ifdef HAVE_OPENMP_OFFLOAD
805 : :
806 : : ABI_MALLOC(work, (mbandc,mbandc*oper%ndat))
807 : : !$OMP TARGET ENTER DATA MAP(alloc:work)
808 : : !do idat=1,oper%ndat,32
809 : : !blk=min(32,oper%ndat-idat+1)
810 : : !$OMP TARGET DATA USE_DEVICE_ADDR(ks,work)
811 : : !call gpu_xginv_strided(2,mbandc,ks(:,1+(idat-1)*mbandc:idat*mbandc,ikpt,isppol),mbandc,mbandc*mbandc,1)
812 : : !call gpu_xginv_strided(2,mbandc,ks(:,1+(idat-1)*mbandc:(idat+blk-1)*mbandc,ikpt,isppol),mbandc,mbandc*mbandc,blk,work)
813 : : call gpu_xginv_strided(2,mbandc,ks(:,:,ikpt,isppol),mbandc,mbandc*mbandc,oper%ndat,work)
814 : : !$OMP END TARGET DATA
815 : : !end do ! idat
816 : : !$OMP TARGET EXIT DATA MAP(delete:work)
817 : : ABI_FREE(work)
818 : :
819 : : #endif
820 : : end if
821 : : end do ! ikpt
822 : : end do ! isppol
823 : : end if ! option
824 : : #ifdef HAVE_OPENMP_OFFLOAD
825 : : !$OMP TARGET EXIT DATA MAP(from:ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
826 : : #endif
827 : :
828 : : ABI_NVTX_END_RANGE()
829 : : DBG_EXIT("COLL")
830 : :
831 [ - + ]: 982154 : end subroutine inverse_oper
832 : : !!***
833 : :
834 : : !!****f* m_oper/downfold_oper
835 : : !! NAME
836 : : !! downfold_oper
837 : : !!
838 : : !! FUNCTION
839 : : !! Downfold an operator from KS space to local space.
840 : : !!
841 : : !! INPUTS
842 : : !! oper <type(oper_type)>= operator
843 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
844 : : !! procb(ikpt)=for kpt parallelization; gives the rank (in the kpt communicator) of the CPU handling each ikpt
845 : : !! iproc=rank of the current process in the kpt communicator
846 : : !! option = 1 (default) : downfold an operator represented by a matrix in KS space
847 : : !! = 2 : downfold the identity
848 : : !! = 3 : downfold a diagonal KS operator
849 : : !! = 4 : computes downfold(upfold)
850 : : !! op_ks_diag = when option=3, you can provide the diagonal KS operator in this variable
851 : : !! with the format mbandc*nkpt*nsppol, instead of storing it in the mband*mband matrix
852 : : !! of oper%ks
853 : : !!
854 : : !! OUTPUT
855 : : !!
856 : : !! SOURCE
857 : :
858 : 302 : subroutine downfold_oper(oper,paw_dmft,procb,iproc,option,op_ks_diag,gpu_option)
859 : :
860 : : !Arguments ------------------------------------
861 : : type(oper_type),target,intent(inout) :: oper
862 : : type(paw_dmft_type),target,intent(in) :: paw_dmft
863 : : integer, optional, intent(in) :: iproc,option,gpu_option
864 : : integer, optional, intent(in) :: procb(oper%nkpt)
865 : : real(dp), optional, intent(in) :: op_ks_diag(oper%mbandc,oper%nkpt,oper%nsppol)
866 : : !oper variables-------------------------------
867 : : integer :: iatom,ib,ik,ikpt,isppol,im,idat,lpawu,mbandc,ndim
868 : : integer :: ndim_max,nspinor,ndat,opt,paral,shift
869 : : integer :: l_gpu_option
870 : : complex(dp) :: alpha
871 : 973940 : complex(dp), ABI_CONTIGUOUS pointer :: ks(:,:,:,:),mat(:,:,:),chipsi(:,:,:,:,:)
872 : 973940 : real(dp), ABI_CONTIGUOUS pointer :: wtk(:)
873 : : character(len=500) :: message
874 : 973940 : complex(dp), target, allocatable :: mat_temp(:,:,:),mat_temp2(:,:,:),mat_temp3(:,:)
875 : : ! *********************************************************************
876 : :
877 : : DBG_ENTER("COLL")
878 : : ABI_NVTX_START_RANGE(NVTX_DMFT_DOWNFOLD_OPER)
879 : :
880 : : #ifndef HAVE_OPENMP_OFFLOAD
881 : : ABI_UNUSED(alpha); ABI_UNUSED(im)
882 : : #endif
883 : :
884 : 973940 : if (oper%has_opermatlu == 0) then
885 : 0 : message = " Operator is not defined to be used in downfold_oper"
886 : 0 : ABI_ERROR(message)
887 : : end if
888 : :
889 [ - + ]: 973940 : l_gpu_option = ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option = gpu_option
890 [ + + + + ]: 973940 : paral = 0; if (present(procb) .and. present(iproc) .and. oper%paral == 0) paral = 1
891 [ + + ]: 973940 : opt = 1; if (present(option)) opt = option
892 : :
893 [ - + ]: 973940 : if(l_gpu_option==ABI_GPU_OPENMP) then
894 [ # # ]: 0 : ABI_CHECK(opt==1 .or. opt==3, "Incompatible codepath with OpenMP GPU")
895 : : end if
896 : :
897 : 973940 : mbandc = oper%mbandc
898 : 973940 : nspinor = oper%nspinor
899 : 973940 : ndim_max = nspinor * (2*paw_dmft%maxlpawu+1)
900 : 973940 : shift = oper%shiftk
901 : 973940 : ndat = oper%ndat
902 : 973940 : ks => oper%ks
903 : 973940 : wtk => oper%wtk
904 : 973940 : chipsi => paw_dmft%chipsi
905 : : if(l_gpu_option==ABI_GPU_OPENMP) then
906 : : #ifdef HAVE_OPENMP_OFFLOAD
907 : : !$OMP TARGET ENTER DATA MAP(to:chipsi,wtk)
908 : : !$OMP TARGET ENTER DATA MAP(to:ks) IF(oper%gpu_option/=ABI_GPU_OPENMP)
909 : : if (present(op_ks_diag)) then
910 : : !$OMP TARGET ENTER DATA MAP(to:op_ks_diag)
911 : : end if
912 : : #endif
913 : : end if
914 : :
915 [ + + ]: 3085141 : do iatom=1,oper%natom
916 : 2111201 : lpawu = oper%matlu(iatom)%lpawu
917 [ + + ]: 2111201 : if (lpawu == -1) cycle
918 : 993976 : mat => oper%matlu(iatom)%mat
919 : 993976 : ndim = nspinor * (2*lpawu+1)
920 [ + - ]: 993976 : if(oper%gpu_option==ABI_GPU_DISABLED) then
921 [ + + + + : 50419989 : mat(:,:,:) = czero
+ + ]
922 [ # # ]: 0 : else if(oper%gpu_option==ABI_GPU_OPENMP) then
923 : 0 : call gpu_set_to_zero_complex(mat, int(oper%nsppol,c_size_t)*ndat*ndim*ndim)
924 : : end if
925 [ - + + - : 4969880 : ABI_MALLOC(mat_temp,(ndim,mbandc,ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
926 [ - + + - : 4969880 : ABI_MALLOC(mat_temp2,(ndim,ndim,ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
927 [ - + + - : 3975904 : ABI_MALLOC(mat_temp3,(ndim,ndim))
- + + - +
- + - - +
- + ]
928 : : #ifdef HAVE_OPENMP_OFFLOAD
929 : : !$OMP TARGET ENTER DATA MAP(alloc:mat_temp,mat_temp2,mat_temp3) IF(l_gpu_option==ABI_GPU_OPENMP)
930 : : !$OMP TARGET ENTER DATA MAP(to:mat) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
931 : : #endif
932 [ + + ]: 2672207 : do isppol=1,oper%nsppol
933 [ + + ]: 11433891 : do ikpt=1,oper%nkpt ! index of kpt on the current CPU
934 : :
935 [ + + ]: 8761684 : if (paral == 1) then
936 [ + + ]: 148596 : if (procb(ikpt) /= iproc) cycle
937 : : end if
938 : :
939 : 8653486 : ik = ikpt + shift ! true kpt index (needed for chipsi)
940 : :
941 [ + + ]: 8653486 : if (opt == 1 .or. opt == 3) then
942 : :
943 : :
944 [ + + ]: 8651270 : if (opt == 1) then
945 : :
946 : : call abi_zgemm_2dd("n","n",ndim,mbandc*ndat,mbandc,cone,&
947 : : & chipsi(:,:,ik,isppol,iatom),ndim_max,&
948 : : & ks(:,:,ikpt,isppol),mbandc,czero,&
949 : : & mat_temp(:,:,:),ndim,&
950 : 8424468 : & gpu_option=l_gpu_option)
951 : :
952 [ + - ]: 226802 : else if (opt == 3) then
953 : :
954 [ + - ]: 226802 : if(l_gpu_option == ABI_GPU_DISABLED) then
955 [ + + ]: 453604 : do idat=1,ndat
956 [ + + ]: 1930938 : do ib=1,mbandc
957 [ + + ]: 1704136 : if (present(op_ks_diag)) then
958 [ + + ]: 109860 : mat_temp(:,ib,idat) = paw_dmft%chipsi(1:ndim,ib,ik,isppol,iatom) * op_ks_diag(ib,ikpt,isppol)
959 : : else
960 [ + + ]: 8852736 : mat_temp(:,ib,idat) = paw_dmft%chipsi(1:ndim,ib,ik,isppol,iatom) * oper%ks(ib,ib+(idat-1)*mbandc,ikpt,isppol)
961 : : end if ! present(op_ks_diag)
962 : : end do ! ib
963 : : end do ! ndat
964 : : else if(l_gpu_option == ABI_GPU_OPENMP) then
965 : : #ifdef HAVE_OPENMP_OFFLOAD
966 : : if (present(op_ks_diag)) then
967 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:chipsi,op_ks_diag,mat_temp) PRIVATE(idat)
968 : : do idat=1,ndat
969 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ib,im)
970 : : do ib=1,mbandc
971 : : do im=1,ndim
972 : : mat_temp(im,ib,idat) = chipsi(im,ib,ik,isppol,iatom) * op_ks_diag(ib,ikpt,isppol)
973 : : end do
974 : : end do
975 : : end do
976 : : else
977 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:chipsi,ks,mat_temp) PRIVATE(idat)
978 : : do idat=1,ndat
979 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ib,im)
980 : : do ib=1,mbandc
981 : : do im=1,ndim
982 : : mat_temp(im,ib,idat) = chipsi(im,ib,ik,isppol,iatom) * ks(ib,ib+(idat-1)*mbandc,ikpt,isppol)
983 : : end do
984 : : end do
985 : : end do
986 : : end if ! present(op_ks_diag)
987 : : #endif
988 : : end if
989 : :
990 : : end if ! opt=1 or 3
991 : :
992 [ + - ]: 8651270 : if(l_gpu_option == ABI_GPU_DISABLED) then
993 [ + + ]: 17302540 : do idat=1,ndat
994 : : call abi_xgemm("n","c",ndim,ndim,mbandc,cone,mat_temp(:,:,idat),ndim,&
995 : 17302540 : & paw_dmft%chipsi(:,:,ik,isppol,iatom),ndim_max,czero,mat_temp2(:,:,idat),ndim)
996 : : end do ! ndat
997 : : else if(l_gpu_option == ABI_GPU_OPENMP) then
998 : : #ifdef HAVE_OPENMP_OFFLOAD
999 : : !$OMP TARGET DATA USE_DEVICE_ADDR(mat_temp,chipsi,mat_temp2)
1000 : : call abi_gpu_xgemm_strided(2,'n','c',ndim,ndim,mbandc,cone,c_loc(mat_temp(:,:,:)),ndim,ndim*mbandc,&
1001 : : & c_loc(chipsi(:,:,ik,isppol,iatom)),ndim_max,0,czero,c_loc(mat_temp2(:,:,:)),ndim,ndim*ndim,ndat)
1002 : : !$OMP END TARGET DATA
1003 : : #endif
1004 : : end if
1005 : :
1006 [ + + ]: 2216 : else if (opt == 2) then
1007 : :
1008 : : call abi_xgemm("n","c",ndim,ndim,mbandc,cone,paw_dmft%chipsi(:,:,ik,isppol,iatom),&
1009 : 1718 : & ndim_max,paw_dmft%chipsi(:,:,ik,isppol,iatom),ndim_max,czero,mat_temp2(:,:,1),ndim)
1010 : :
1011 [ + - ]: 498 : else if (opt == 4) then
1012 : :
1013 : : call abi_xgemm("n","c",ndim,ndim,mbandc,cone,paw_dmft%chipsi(:,:,ik,isppol,iatom),&
1014 : 498 : & ndim_max,paw_dmft%chipsi(:,:,ik,isppol,iatom),ndim_max,czero,mat_temp3(:,:),ndim)
1015 : :
1016 : : call abi_xgemm("n","n",ndim,ndim,ndim,cone,mat_temp3(:,:),ndim,&
1017 : 498 : & mat_temp3(:,:),ndim,czero,mat_temp2(:,:,1),ndim)
1018 : :
1019 : : end if ! opt
1020 : :
1021 [ + - ]: 10331717 : if(l_gpu_option == ABI_GPU_DISABLED) then
1022 [ + + ]: 17306972 : do idat=1,ndat
1023 : : oper%matlu(iatom)%mat(:,:,idat+(isppol-1)*ndat) = &
1024 [ + + + + ]: 272685924 : & oper%matlu(iatom)%mat(:,:,idat+(isppol-1)*ndat) + mat_temp2(:,:,idat)*oper%wtk(ik)
1025 : : end do ! ndat
1026 : : else if(l_gpu_option == ABI_GPU_OPENMP) then
1027 : : #ifdef HAVE_OPENMP_OFFLOAD
1028 : : alpha = dcmplx(wtk(ik), 0.0_dp)
1029 : : call abi_xaxpy(ndim*ndim*ndat, alpha, &
1030 : : & mat_temp2, 1, mat(:,:,1+(isppol-1)*ndat:isppol*ndat), 1, gpu_option=l_gpu_option)
1031 : : #endif
1032 : : end if
1033 : :
1034 : : end do ! ikpt
1035 : : end do ! isppol
1036 : : #ifdef HAVE_OPENMP_OFFLOAD
1037 : : !$OMP TARGET EXIT DATA MAP(from:mat) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1038 : : !$OMP TARGET EXIT DATA MAP(delete:mat_temp,mat_temp2,mat_temp3) IF(l_gpu_option==ABI_GPU_OPENMP)
1039 : : #endif
1040 [ - + ]: 993976 : ABI_FREE(mat_temp)
1041 [ - + ]: 993976 : ABI_FREE(mat_temp2)
1042 [ - + ]: 3085141 : ABI_FREE(mat_temp3)
1043 : : end do ! iatom
1044 : :
1045 : : #ifdef HAVE_OPENMP_OFFLOAD
1046 : : !$OMP TARGET EXIT DATA MAP(delete:ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1047 : : !$OMP TARGET EXIT DATA MAP(delete:chipsi,wtk) IF(l_gpu_option==ABI_GPU_OPENMP)
1048 : : if (present(op_ks_diag)) then
1049 : : !$OMP TARGET EXIT DATA MAP(delete:op_ks_diag) IF(l_gpu_option==ABI_GPU_OPENMP)
1050 : : end if
1051 : : #endif
1052 : : !do isppol=1,nsppol
1053 : : ! do ikpt=1,nkpt
1054 : : ! ikpt1=ikpt
1055 : : ! if(present(jkpt)) ikpt1=jkpt
1056 : : ! lvz=paral==0 !vz_d
1057 : : ! if(present(iproc)) lvz=lvz.or.(paral==1.and.(procb2(ikpt1)==iproc)) !vz_d
1058 : : !! if ((paral==1.and.(procb2(ikpt1)==iproc)).or.(paral==0)) then !vz_d
1059 : : ! if(lvz) then !vz_d
1060 : : ! do ib=1,mbandc
1061 : : ! do ib1=1,mbandc
1062 : : ! do iatom=1,natom
1063 : : ! if(oper%matlu(iatom)%lpawu.ne.-1) then
1064 : : ! ndim=2*oper%matlu(iatom)%lpawu+1
1065 : : ! do im=1,ndim
1066 : : ! do im1=1,ndim
1067 : : ! do ispinor=1,nspinor
1068 : : ! do ispinor1=1,nspinor
1069 : : ! if (im1 == im .and. im1 == 1 .and. ib1 == ib .and. iatom == 1) then
1070 : :
1071 : : ! end if
1072 : : ! oper%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)= &
1073 : : !& oper%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)+ &
1074 : : !& paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)* &
1075 : : !& conjg(paw_dmft%psichi(isppol,ikpt1,ib1,ispinor1,iatom,im1))* &
1076 : : !false& paw_dmft%psichi(isppol,ikpt1,ib1,ispinor1,iatom,im1)*
1077 : : !&
1078 : : !false& conjg(paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im))* &
1079 : : !& oper%ks(isppol,ikpt,ib,ib1)*oper%wtk(ikpt)
1080 : : ! one could suppress wtk here if present(jkpt)
1081 : : ! ks(ib,ib1)=ks(ib1,ib) -> ib and ib1 can be underchanged !
1082 : : ! enddo ! ispinor1
1083 : : ! enddo ! ispinor
1084 : : ! enddo ! im1
1085 : : ! enddo ! im
1086 : : ! endif
1087 : : ! enddo ! iatom
1088 : : ! enddo ! ib
1089 : : ! enddo ! ib
1090 : : ! endif
1091 : : ! enddo ! ikpt
1092 : : !enddo ! isppol
1093 : :
1094 : : DBG_EXIT("COLL")
1095 : :
1096 : : ABI_NVTX_END_RANGE()
1097 [ + + - + ]: 974242 : end subroutine downfold_oper
1098 : : !!***
1099 : :
1100 : :
1101 : : !!****f* m_oper/upfold_oper
1102 : : !! NAME
1103 : : !! upfold_oper
1104 : : !!
1105 : : !! FUNCTION
1106 : : !! Upfold an operator from local space to KS space
1107 : : !!
1108 : : !! INPUTS
1109 : : !! oper <type(oper_type)>= operator
1110 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1111 : : !! procb(ikpt)=for kpt parallelization; gives the rank (in the kpt communicator) of the CPU handling each ikpt
1112 : : !! iproc=rank of the current process in the kpt communicator
1113 : : !!
1114 : : !! OUTPUT
1115 : : !!
1116 : : !! SOURCE
1117 : :
1118 : 949905 : subroutine upfold_oper(oper,paw_dmft,procb,iproc,gpu_option)
1119 : :
1120 : : !Arguments ------------------------------------
1121 : : type(oper_type),target, intent(inout) :: oper
1122 : : type(paw_dmft_type),target, intent(in) :: paw_dmft
1123 : : integer, optional, intent(in) :: iproc,gpu_option
1124 : : integer, optional, intent(in) :: procb(oper%nkpt)
1125 : : !Local variables-------------------------------
1126 : : integer :: iatom,ik,ikpt,isppol,idat,lpawu,mbandc,l_gpu_option
1127 : : integer :: ndim,ndim_max,ndat,nspinor,paral,shift
1128 : 949905 : complex(dp), ABI_CONTIGUOUS pointer :: ks(:,:,:,:),mat(:,:,:),chipsi(:,:,:,:,:)
1129 : 949905 : complex(dp), target, allocatable :: mat_temp(:,:),mat_temp2(:,:)
1130 : : ! *********************************************************************
1131 : :
1132 : 0 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option=gpu_option
1133 : :
1134 : : DBG_ENTER("COLL")
1135 : : ABI_NVTX_START_RANGE(NVTX_DMFT_UPFOLD_OPER)
1136 : :
1137 : : !if ((oper%has_opermatlu == 0) .or. (oper%has_operks == 0)) then
1138 : : ! message = " Operator is not defined to be used in upfold_oper"
1139 : : ! ABI_ERROR(message)
1140 : : !end if
1141 : :
1142 : 949905 : mbandc = paw_dmft%mbandc
1143 : 949905 : nspinor = paw_dmft%nspinor
1144 : 949905 : ndim_max = nspinor * (2*paw_dmft%maxlpawu+1)
1145 : 949905 : paral = 0
1146 : 949905 : shift = oper%shiftk
1147 : 949905 : ndat = oper%ndat
1148 : :
1149 [ - + + + ]: 949905 : if (present(procb) .and. present(iproc) .and. oper%paral == 0) paral = 1
1150 : :
1151 [ + - + + : 504189235 : if(l_gpu_option==ABI_GPU_DISABLED) oper%ks(:,:,:,:) = czero
+ + + + +
+ ]
1152 : :
1153 [ - + + - : 3799620 : ABI_MALLOC(mat_temp,(mbandc,ndim_max*ndat))
- + + - +
- + - - +
- + ]
1154 [ - + + - : 3799620 : ABI_MALLOC(mat_temp2,(mbandc,mbandc*ndat))
- + + - +
- + - - +
- + ]
1155 : 949905 : ks => oper%ks
1156 : 949905 : chipsi => paw_dmft%chipsi
1157 : :
1158 : : if(l_gpu_option == ABI_GPU_OPENMP) then
1159 : : #ifdef HAVE_OPENMP_OFFLOAD
1160 : : !$OMP TARGET ENTER DATA MAP(alloc:ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1161 : : !$OMP TARGET ENTER DATA MAP(alloc:chipsi,mat_temp,mat_temp2) IF(l_gpu_option==ABI_GPU_OPENMP)
1162 : : !$OMP TARGET UPDATE TO(chipsi) IF(l_gpu_option==ABI_GPU_OPENMP)
1163 : : call gpu_set_to_zero_complex(ks, int(oper%nsppol,c_size_t)*ndat*mbandc*mbandc*oper%nkpt)
1164 : : #endif
1165 : : end if
1166 : :
1167 [ + + ]: 3014050 : do iatom=1,oper%natom
1168 : 2064145 : lpawu = oper%matlu(iatom)%lpawu
1169 [ + + ]: 2064145 : if (lpawu == -1) cycle
1170 : 969393 : ndim = (2*lpawu+1) * nspinor
1171 : 969393 : mat => oper%matlu(iatom)%mat
1172 : : #ifdef HAVE_OPENMP_OFFLOAD
1173 : : !$OMP TARGET ENTER DATA MAP(to:mat) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1174 : : #endif
1175 [ + + ]: 3554624 : do isppol=1,oper%nsppol
1176 [ + + ]: 12218723 : do ikpt=1,oper%nkpt ! index of kpt on the current CPU
1177 : :
1178 [ + + ]: 8519252 : if (paral == 1) then
1179 [ + + ]: 132000 : if (procb(ikpt) /= iproc) cycle
1180 : : end if
1181 : :
1182 : 8420252 : ik = ikpt + shift ! true kpt index (needed for chipsi)
1183 : :
1184 [ + - ]: 10055578 : if(l_gpu_option == ABI_GPU_DISABLED) then
1185 : :
1186 : : call abi_zgemm_2dd("c","n",mbandc,ndat*ndim,ndim,cone,paw_dmft%chipsi(:,:,ik,isppol,iatom),&
1187 : 8420252 : & ndim_max,oper%matlu(iatom)%mat(:,:,(isppol-1)*ndat+1:isppol*ndat),ndim,czero,mat_temp(:,:),mbandc)
1188 : :
1189 [ + + ]: 16840504 : do idat=1,ndat
1190 : :
1191 : : call abi_xgemm("n","n",mbandc,mbandc,ndim,cone,mat_temp(:,1+(idat-1)*ndim:idat*ndim),mbandc,&
1192 : 16840504 : & paw_dmft%chipsi(:,:,ik,isppol,iatom),ndim_max,czero,mat_temp2(:,1+(idat-1)*mbandc:idat*mbandc),mbandc)
1193 : :
1194 : : end do ! idat
1195 : :
1196 : : !oper%ks(:,:,ikpt,isppol) = oper%ks(:,:,ikpt,isppol) + mat_temp2(:,:)
1197 : 8420252 : call zaxpy(mbandc*mbandc*ndat, cone, mat_temp2, 1, oper%ks(:,:,ikpt,isppol), 1)
1198 : :
1199 : : else if(l_gpu_option == ABI_GPU_OPENMP) then
1200 : : #ifdef HAVE_OPENMP_OFFLOAD
1201 : : call abi_zgemm("c","n",mbandc,ndat*ndim,ndim,cone,&
1202 : : & chipsi(:,:,ik,isppol,iatom),ndim_max,&
1203 : : & mat(:,:,(isppol-1)*ndat+1:isppol*ndat),ndim,czero,&
1204 : : & mat_temp(:,:),mbandc,&
1205 : : & gpu_option=l_gpu_option)
1206 : :
1207 : : !$OMP TARGET DATA USE_DEVICE_ADDR(mat_temp,chipsi,mat_temp2)
1208 : : call abi_gpu_xgemm_strided(2,'n','n',mbandc,mbandc,ndim,cone,c_loc(mat_temp(:,:)),mbandc,ndim*mbandc,&
1209 : : & c_loc(chipsi(:,:,ik,isppol,iatom)),ndim_max,0,czero,c_loc(mat_temp2(:,:)),mbandc,mbandc*mbandc,ndat)
1210 : : !$OMP END TARGET DATA
1211 : : !$OMP TARGET DATA USE_DEVICE_ADDR(ks,mat_temp2)
1212 : : call abi_gpu_xaxpy(1, 2*mbandc*mbandc*ndat, cone, &
1213 : : & c_loc(mat_temp2), 1, c_loc(ks(:,:,ikpt,isppol)), 1)
1214 : : !$OMP END TARGET DATA
1215 : : #endif
1216 : : end if
1217 : :
1218 : : end do ! ikpt
1219 : : end do ! isppol
1220 : : #ifdef HAVE_OPENMP_OFFLOAD
1221 : : !$OMP TARGET EXIT DATA MAP(delete:mat) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1222 : : #endif
1223 : : end do ! iatom
1224 : :
1225 : : #ifdef HAVE_OPENMP_OFFLOAD
1226 : : !$OMP TARGET UPDATE FROM(ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1227 : : !$OMP TARGET EXIT DATA MAP(delete:ks) IF(l_gpu_option==ABI_GPU_OPENMP .and. oper%gpu_option/=ABI_GPU_OPENMP)
1228 : : !$OMP TARGET EXIT DATA MAP(delete:chipsi,mat_temp,mat_temp2) IF(l_gpu_option==ABI_GPU_OPENMP)
1229 : : #endif
1230 [ - + ]: 949905 : ABI_FREE(mat_temp)
1231 [ - + ]: 949905 : ABI_FREE(mat_temp2)
1232 : :
1233 : : ABI_NVTX_END_RANGE()
1234 : : DBG_EXIT("COLL")
1235 : :
1236 [ - + ]: 949905 : end subroutine upfold_oper
1237 : : !!***
1238 : :
1239 : : !!****f* m_oper/identity_oper
1240 : : !! NAME
1241 : : !! identity_oper
1242 : : !!
1243 : : !! FUNCTION
1244 : : !! Construct the identity operator
1245 : : !!
1246 : : !! INPUTS
1247 : : !! oper <type(oper_type)>= operator
1248 : : !! option = 1: in KS space
1249 : : !! = 2: in local space
1250 : : !! = 3: both
1251 : : !!
1252 : : !! OUTPUT
1253 : : !!
1254 : : !! SOURCE
1255 : :
1256 : 175 : subroutine identity_oper(oper,option)
1257 : :
1258 : : !Arguments ------------------------------------
1259 : : integer, intent(in) :: option
1260 : : type(oper_type), intent(inout) :: oper
1261 : : !Local variables-------------------------------
1262 : : integer :: ib,natom
1263 : : character(len=500) :: message
1264 : : ! *********************************************************************
1265 : :
1266 : : DBG_ENTER("COLL")
1267 : :
1268 [ - + - - : 175 : if (((option == 1 .or. option == 3) .and. (oper%has_operks == 0)) .or. &
+ - - + ]
1269 : : & ((option == 2 .or. option == 3) .and. (oper%has_opermatlu == 0))) then
1270 : 0 : message = " Options in identity_oper are not coherent with definitions of this operator"
1271 : 0 : ABI_ERROR(message)
1272 : : end if
1273 : :
1274 [ - + ]: 175 : if (option == 1 .or. option == 3) then
1275 : :
1276 [ # # # # : 0 : oper%ks(:,:,:,:) = czero
# # # # ]
1277 [ # # ]: 0 : do ib=1,oper%mbandc
1278 [ # # # # ]: 0 : oper%ks(ib,ib,:,:) = cone
1279 : : end do ! ib
1280 : :
1281 : : end if ! option=1 or 3
1282 : :
1283 [ + - ]: 175 : if (option == 2 .or. option == 3) then
1284 : 175 : natom = oper%natom
1285 : 175 : call zero_matlu(oper%matlu(:),natom)
1286 : 175 : call identity_matlu(oper%matlu(:),natom)
1287 : : end if ! option=2 or 3
1288 : :
1289 : : DBG_EXIT("COLL")
1290 : :
1291 : 175 : end subroutine identity_oper
1292 : : !!***
1293 : :
1294 : : !!****f* m_oper/diff_oper
1295 : : !! NAME
1296 : : !! diff_oper
1297 : : !!
1298 : : !! FUNCTION
1299 : : !! Compute a norm of the differences between two occupations matrices.
1300 : : !!
1301 : : !! INPUTS
1302 : : !! char1 = character describing occup1
1303 : : !! char2 = character describing occup2
1304 : : !! occup1 <type(oper_type)>= occupations
1305 : : !! occup2 <type(oper_type)>= occupations
1306 : : !! option : option for printing (if 1 assume data are related to lda only)
1307 : : !! toldiff : tolerance for the difference
1308 : : !!
1309 : : !! OUTPUT
1310 : : !!
1311 : : !! SOURCE
1312 : :
1313 : 172 : subroutine diff_oper(char1,char2,occup1,occup2,option,toldiff)
1314 : :
1315 : : !Arguments ------------------------------------
1316 : : type(oper_type), intent(in) :: occup1,occup2
1317 : : integer, intent(in) :: option
1318 : : real(dp), intent(in) :: toldiff
1319 : : character(len=*), intent(in) :: char1,char2
1320 : : !Local variables-------------------------------
1321 : : character(len=500) :: message
1322 : : ! *********************************************************************
1323 : :
1324 : : DBG_ENTER("COLL")
1325 : :
1326 [ + - - + ]: 172 : if (occup1%has_opermatlu == 0 .or. occup2%has_opermatlu == 0) then
1327 : 0 : message = " Operators are not defined to be used in diff_oper"
1328 : 0 : ABI_ERROR(message)
1329 : : end if
1330 : :
1331 [ - + ]: 172 : if (occup1%nkpt /= occup2%nkpt) then
1332 : 0 : write(message,'(a,2x,2i9)')' Operators are not equal',occup1%nkpt,occup2%nkpt
1333 : 0 : ABI_ERROR(message)
1334 : : end if
1335 : :
1336 : 172 : call diff_matlu(char1,char2,occup1%matlu(:),occup2%matlu(:),occup1%natom,option,toldiff)
1337 : : ! if(option==1) then
1338 : : ! toldiff=tol4
1339 : : ! if( matludiff < toldiff ) then
1340 : : ! write(message,'(6a,e12.4,a,e12.4)') ch10,&
1341 : : !& ' Differences between ',trim(char1),' and ',trim(char2),' is small enough:',&
1342 : : !& matludiff,'is lower than',toldiff
1343 : : ! call wrtout(std_out,message,'COLL')
1344 : : ! else
1345 : : ! write(message,'(6a,e12.4,a,e12.4)') ch10,&
1346 : : !& ' Error: Differences between ',trim(char1),' and ',trim(char2),' is too large:',&
1347 : : !& matludiff,'is large than',toldiff
1348 : : ! call wrtout(std_out,message,'COLL')
1349 : : ! call abi_abort('COLL')
1350 : : ! endif
1351 : : ! endif
1352 : : ! call abi_abort('COLL')
1353 : :
1354 : : DBG_EXIT("COLL")
1355 : :
1356 : 172 : end subroutine diff_oper
1357 : : !!***
1358 : :
1359 : : !!****f* m_oper/trace_oper
1360 : : !! NAME
1361 : : !! trace_oper
1362 : : !!
1363 : : !! FUNCTION
1364 : : !! Computes the trace of an operator
1365 : : !!
1366 : : !! INPUTS
1367 : : !! oper <type(oper_type)>= operator
1368 : : !! opt_ksloc = 1: trace in KS space
1369 : : !! = 2: trace in local space
1370 : : !! = 3: both
1371 : : !!
1372 : : !! OUTPUT
1373 : : !! trace_ks :: trace in KS space
1374 : : !! trace_loc :: trace in local space
1375 : : !! trace_ks_cmplx :: complex trace in KS space
1376 : : !!
1377 : : !! SOURCE
1378 : :
1379 : 5659 : subroutine trace_oper(oper,trace_ks,trace_loc,opt_ksloc,trace_ks_cmplx)
1380 : :
1381 : : !Arguments ------------------------------------
1382 : : type(oper_type), intent(in) :: oper
1383 : : real(dp), intent(out) :: trace_ks !vz_i
1384 : : real(dp), intent(inout) :: trace_loc(oper%nsppol+1,oper%natom) !vz_i
1385 : : integer, intent(in) :: opt_ksloc
1386 : : complex(dp), optional, intent(out) :: trace_ks_cmplx
1387 : : !Local variables-------------------------------
1388 : : integer :: ib,ikpt,isppol
1389 : : complex(dp) :: trace
1390 : : character(len=500) :: message
1391 : : ! *********************************************************************
1392 : :
1393 : : DBG_ENTER("COLL")
1394 : :
1395 [ + + + - : 5659 : if (((opt_ksloc == 1 .or. opt_ksloc == 3) .and. (oper%has_operks == 0)) .or. &
+ + - + ]
1396 : : & ((opt_ksloc == 2 .or. opt_ksloc == 3) .and. (oper%has_opermatlu == 0))) then
1397 : 0 : message = " Options in trace_oper are not coherent with definitions of this operator"
1398 : 0 : ABI_ERROR(message)
1399 : : end if
1400 : :
1401 [ + + ]: 5659 : if (opt_ksloc == 1 .or. opt_ksloc == 3) then
1402 : 4880 : trace = czero
1403 : : !temp1=zero
1404 [ + + ]: 11564 : do isppol=1,oper%nsppol
1405 [ + + ]: 60270 : do ikpt=1,oper%nkpt
1406 [ + + ]: 534886 : do ib=1,oper%mbandc
1407 : 528202 : trace = trace + oper%ks(ib,ib,ikpt,isppol)*oper%wtk(ikpt+oper%shiftk)
1408 : : !temp1=temp1+oper%wtk(ikpt)
1409 : : end do ! ib
1410 : : end do ! ikpt
1411 : : end do ! isppol
1412 [ + + + + ]: 4880 : if (oper%nsppol == 1 .and. oper%nspinor == 1) trace = two * trace
1413 [ - + ]: 4880 : if (present(trace_ks_cmplx)) trace_ks_cmplx = trace
1414 : 4880 : trace_ks = dble(trace)
1415 : : ! write(std_out,*) "temp1",temp1
1416 : : end if ! opt_ksloc
1417 : :
1418 [ + + ]: 5659 : if (opt_ksloc == 2 .or. opt_ksloc == 3) then
1419 : 779 : call trace_matlu(oper%matlu(:),oper%natom,trace_loc=trace_loc(:,:))
1420 : : end if
1421 : :
1422 : : DBG_EXIT("COLL")
1423 : :
1424 : 5659 : end subroutine trace_oper
1425 : : !!***
1426 : :
1427 : : !!****f* m_oper/prod_oper
1428 : : !! NAME
1429 : : !! prod_oper
1430 : : !!
1431 : : !! FUNCTION
1432 : : !! Computes the matrix product of oper1 and oper2
1433 : : !!
1434 : : !! INPUTS
1435 : : !! oper1,oper2 <type(oper_type)>= operator
1436 : : !! opt_ksloc = 1 : in KS space
1437 : : !! = 2 : in local space
1438 : : !! opt_diag = 1 if oper1 and oper2 are diagonal in KS space, 0 otherwise (default)
1439 : : !!
1440 : : !! OUTPUT
1441 : : !! oper3 <type(oper_type)>= matrix product of oper1 and oper2
1442 : : !!
1443 : : !! SOURCE
1444 : :
1445 : 68 : subroutine prod_oper(oper1,oper2,oper3,opt_ksloc,opt_diag)
1446 : :
1447 : : !Arguments ------------------------------------
1448 : : type(oper_type), intent(in) :: oper1,oper2
1449 : : type(oper_type), intent(inout) :: oper3
1450 : : integer, intent(in) :: opt_ksloc
1451 : : integer, optional, intent(in) :: opt_diag
1452 : : !Local variables-------------------------------
1453 : : integer :: ib,ikpt,isppol,mbandc
1454 : : logical :: diag
1455 : : ! *********************************************************************
1456 : :
1457 : : DBG_ENTER("COLL")
1458 : :
1459 : : if (opt_ksloc == 2 .and. oper1%has_opermatlu == 1 .and. &
1460 [ + + + - : 68 : & oper2%has_opermatlu == 1 .and. oper3%has_opermatlu == 1) then
+ - + - ]
1461 : 34 : call prod_matlu(oper1%matlu(:),oper2%matlu(:),oper3%matlu(:),oper1%natom)
1462 : : end if
1463 : :
1464 : : if (opt_ksloc == 1 .and. oper1%has_operks == 1 .and. &
1465 [ + + + - : 68 : & oper2%has_operks == 1 .and. oper3%has_operks == 1) then
+ - + - ]
1466 : 34 : mbandc = oper1%mbandc
1467 : 34 : diag = .false.
1468 [ - + ]: 34 : if (present(opt_diag)) then
1469 [ # # ]: 0 : if (opt_diag == 1) diag = .true.
1470 : : end if
1471 : : if (diag) then
1472 [ # # ]: 0 : do ib=1,mbandc
1473 [ # # # # ]: 0 : oper3%ks(ib,ib,:,:) = oper1%ks(ib,ib,:,:) * oper2%ks(ib,ib,:,:)
1474 : : end do ! ib
1475 : : else
1476 [ + + ]: 88 : do isppol=1,oper1%nsppol
1477 [ + + ]: 888 : do ikpt=1,oper1%nkpt
1478 : : call abi_xgemm("n","n",mbandc,mbandc,mbandc,cone,oper1%ks(:,:,ikpt,isppol),mbandc,&
1479 : 854 : & oper2%ks(:,:,ikpt,isppol),mbandc,czero,oper3%ks(:,:,ikpt,isppol),mbandc)
1480 : : end do ! ikpt
1481 : : end do ! isppol
1482 : : end if ! diag
1483 : : end if ! opt_ksloc=1
1484 : :
1485 : : DBG_EXIT("COLL")
1486 : :
1487 : 68 : end subroutine prod_oper
1488 : : !!***
1489 : :
1490 : : !!****f* m_oper/trace_prod_oper
1491 : : !! NAME
1492 : : !! trace_prod_oper
1493 : : !!
1494 : : !! FUNCTION
1495 : : !! Computes Tr(oper1*oper2) in KS space
1496 : : !!
1497 : : !! INPUTS
1498 : : !! oper1,oper2 <type(oper_type)>= operator
1499 : : !!
1500 : : !! OUTPUT
1501 : : !! trace = Tr(op1*op2)
1502 : : !!
1503 : : !! SOURCE
1504 : :
1505 : 0 : subroutine trace_prod_oper(oper1,oper2,trace)
1506 : :
1507 : : !Arguments ------------------------------------
1508 : : type(oper_type), intent(in) :: oper1,oper2
1509 : : complex(dp), intent(out) :: trace
1510 : : !Local variables-------------------------------
1511 : : integer :: ikpt,isppol
1512 : : character(len=500) :: message
1513 : : ! *********************************************************************
1514 : :
1515 [ # # ]: 0 : if (oper1%shiftk /= oper2%shiftk) then
1516 : 0 : message = "Inconsistency in trace_prod_oper: oper1%shiftk should be equal to oper2%shiftk"
1517 : 0 : ABI_ERROR(message)
1518 : : end if
1519 : :
1520 : 0 : trace = czero
1521 : :
1522 [ # # ]: 0 : do isppol=1,oper1%nsppol
1523 [ # # ]: 0 : do ikpt=1,oper1%nkpt
1524 : : trace = trace + sum(oper1%ks(:,:,ikpt,isppol)*transpose(oper2%ks(:,:,ikpt,isppol)))*&
1525 [ # # # # ]: 0 : & oper1%wtk(ikpt+oper1%shiftk)
1526 : : end do ! ikpt
1527 : : end do ! isppol
1528 : :
1529 [ # # # # ]: 0 : if (oper1%nsppol == 1 .and. oper1%nspinor == 1) trace = trace * two
1530 : :
1531 : 0 : end subroutine trace_prod_oper
1532 : : !!***
1533 : :
1534 : : !!****f* m_oper/gather_oper
1535 : : !! NAME
1536 : : !! gather_oper
1537 : : !!
1538 : : !! FUNCTION
1539 : : !! Gather the contributions from all CPUs, for a frequency-dependent
1540 : : !! operator, and for both levels of parallelization (kpt and then frequency,
1541 : : !! and frequency only).
1542 : : !!
1543 : : !! INPUTS
1544 : : !! oper <type(oper_type)>= operator for each frequency
1545 : : !! distrib <type(mpi_distrib_dmft_type)> = mpi related data
1546 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1547 : : !! opt_ksloc = 1 : gather the KS operator on the kpt and frequency communicator
1548 : : !! 2 : gather the local operator (the exact behavior can be defined via opt_commkpt)
1549 : : !! master = if present, only gather on the master node
1550 : : !! opt_diag = 1 if the operator is diagonal in KS space, 0 (default) otherwise
1551 : : !! opt_commkpt (only meaningful for the local quantity)
1552 : : !! = 0 (default) : frequency-only parallelization
1553 : : !! -> xmpi_allgatherv on the whole communicator
1554 : : !! = 1 : kpt and then frequency parallelization (CAREFUL: here, frequencies
1555 : : !! are not distributed in the same way as the frequency-only parallelization scheme)
1556 : : !! -> xmpi_sum on the kpt-communicator, and then xmpi_allgatherv
1557 : : !! on the frequency communicator (useful after a downfold for instance)
1558 : : !!
1559 : : !! OUTPUT
1560 : : !!
1561 : : !! SOURCE
1562 : :
1563 : 1506 : subroutine gather_oper(oper,distrib,paw_dmft,opt_ksloc,master,opt_diag,opt_commkpt)
1564 : :
1565 : : !Arguments ------------------------------------
1566 : : type(mpi_distrib_dmft_type), target, intent(in) :: distrib
1567 : : type(oper_type), intent(inout) :: oper(distrib%nw)
1568 : : type(paw_dmft_type) :: paw_dmft
1569 : : integer, intent(in) :: opt_ksloc
1570 : : integer, optional, intent(in) :: master,opt_commkpt,opt_diag
1571 : : !Local variables-------------------------------
1572 : : integer :: comm,iatom,ib1,ibuf,ierr,ifreq,ikpt,im1
1573 : : integer :: irank,irank1,irank2,isppol,lpawu,mbandc,myproc
1574 : : integer :: myproc2,natom,ndim,nkpt,nproc,nproc_freq,nproc_kpt
1575 : : integer :: nproc2,nspinor,nsppol,nw,optcommkpt,siz_buf
1576 : : logical :: diag
1577 : 1506 : integer, allocatable :: displs(:),recvcounts(:)
1578 : 1506 : complex(dp), allocatable :: buffer(:),buffer_tot(:)
1579 : : ! *********************************************************************
1580 : :
1581 : 1506 : comm = paw_dmft%spacecomm
1582 : 1506 : mbandc = paw_dmft%mbandc
1583 : 1506 : myproc = paw_dmft%myproc
1584 : 1506 : natom = paw_dmft%natom
1585 : 1506 : nkpt = paw_dmft%nkpt
1586 : 1506 : nproc = paw_dmft%nproc
1587 : 1506 : nspinor = paw_dmft%nspinor
1588 : 1506 : nsppol = paw_dmft%nsppol
1589 : 1506 : nw = distrib%nw
1590 : :
1591 : 1506 : nproc_kpt = min(nkpt,nproc)
1592 : :
1593 : 1506 : optcommkpt = 0
1594 [ + + ]: 1506 : if (present(opt_commkpt)) optcommkpt = opt_commkpt
1595 : :
1596 [ + + ]: 1506 : if (opt_ksloc == 1) then
1597 : :
1598 : 4 : nproc_freq = max(1,nproc/nkpt)
1599 : :
1600 [ - + + - : 12 : ABI_MALLOC(recvcounts,(nproc))
+ - + - -
+ - + ]
1601 [ - + + - : 8 : ABI_MALLOC(displs,(nproc))
- + - + ]
1602 : :
1603 : 4 : diag = .false.
1604 [ - + ]: 4 : if (present(opt_diag)) then
1605 [ + - ]: 4 : if (opt_diag == 1) diag = .true.
1606 : : end if
1607 : :
1608 : 4 : irank2 = 1
1609 [ + + ]: 20 : do irank=0,nproc_kpt-1
1610 [ + + ]: 36 : do irank1=0,nproc_freq-1
1611 : 16 : recvcounts(irank2) = distrib%nkpt_mem(irank+1) * distrib%nw_mem_kptparal(irank1+1)
1612 : 32 : irank2 = irank2 + 1
1613 : : end do ! irank1
1614 : : end do ! irank
1615 [ - + - - ]: 4 : if (nproc > nproc_freq*nproc_kpt) recvcounts(nproc_freq*nproc_kpt+1:nproc) = 0
1616 : :
1617 [ + - + + ]: 24 : recvcounts(:) = recvcounts(:) * merge(mbandc,mbandc**2,diag)
1618 : 4 : displs(1) = 0
1619 [ + + ]: 16 : do irank=2,nproc
1620 : 16 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
1621 : : end do ! irank
1622 : :
1623 [ - + + - : 12 : ABI_MALLOC(buffer,(recvcounts(myproc+1)))
+ - + - -
+ - + ]
1624 [ - + + - : 12 : ABI_MALLOC(buffer_tot,(displs(nproc)+recvcounts(nproc)))
+ - + - -
+ - + ]
1625 : :
1626 [ + + ]: 8 : do isppol=1,nsppol
1627 : :
1628 : : ibuf = 0
1629 [ + + ]: 48 : do ikpt=1,nkpt
1630 [ + + ]: 44 : if (distrib%procb(ikpt) /= distrib%me_kpt) cycle
1631 [ + + ]: 33015 : do ifreq=1,nw
1632 [ + - ]: 33000 : if (distrib%proct(ifreq) /= distrib%me_freq) cycle
1633 [ + + ]: 330044 : do ib1=1,mbandc
1634 [ - + ]: 330000 : if (diag) then
1635 : 0 : ibuf = ibuf + 1
1636 : 0 : buffer(ibuf) = oper(ifreq)%ks(ib1,ib1,ikpt,isppol)
1637 : : else
1638 [ + + ]: 2970000 : buffer(ibuf+1:ibuf+mbandc) = oper(ifreq)%ks(:,ib1,ikpt,isppol)
1639 : : ibuf = ibuf + mbandc
1640 : : end if ! diag
1641 : : end do ! ib1
1642 : : end do ! ifreq
1643 : : end do ! ikpt
1644 : :
1645 [ - + ]: 4 : if (present(master)) then
1646 : 0 : call xmpi_gatherv(buffer(:),recvcounts(myproc+1),buffer_tot(:),recvcounts(:),displs(:),master,comm,ierr)
1647 : : else
1648 : 4 : call xmpi_allgatherv(buffer(:),recvcounts(myproc+1),buffer_tot(:),recvcounts(:),displs(:),comm,ierr)
1649 : : end if ! present(master)
1650 : :
1651 : : ibuf = 0
1652 [ + + ]: 52 : do ikpt=1,nkpt
1653 [ + + ]: 132048 : do ifreq=1,nw
1654 [ + + ]: 1320044 : do ib1=1,mbandc
1655 [ - + ]: 1320000 : if (diag) then
1656 : 0 : ibuf = ibuf + 1
1657 : 0 : oper(ifreq)%ks(ib1,ib1,ikpt,isppol) = buffer_tot(ibuf)
1658 : : else
1659 [ + + ]: 11880000 : oper(ifreq)%ks(:,ib1,ikpt,isppol) = buffer_tot(ibuf+1:ibuf+mbandc)
1660 : : ibuf = ibuf + mbandc
1661 : : end if ! diag
1662 : : end do ! ib1
1663 : : end do ! ifreq
1664 : : end do ! ikpt
1665 : :
1666 : : end do ! isppol
1667 : :
1668 [ + - ]: 1502 : else if (opt_ksloc == 2) then
1669 : :
1670 : 1502 : nproc_freq = nproc / nkpt
1671 [ + + ]: 1502 : myproc2 = merge(distrib%me_freq,myproc,optcommkpt==1)
1672 : 587 : nproc2 = merge(nproc_freq+1,nproc,optcommkpt==1)
1673 : :
1674 [ - + + - : 4506 : ABI_MALLOC(recvcounts,(nproc2))
+ - + - -
+ - + ]
1675 [ - + + - : 3004 : ABI_MALLOC(displs,(nproc2))
- + + - ]
1676 : :
1677 : 6672 : siz_buf = 0
1678 : :
1679 [ + + ]: 6672 : do iatom=1,natom
1680 : 5170 : lpawu = paw_dmft%lpawu(iatom)
1681 [ + + ]: 5170 : if (lpawu == -1) cycle
1682 : 6672 : siz_buf = siz_buf + (2*lpawu+1)**2
1683 : : end do ! iatom
1684 : :
1685 : 1502 : siz_buf = siz_buf * (nspinor**2) * nsppol
1686 [ + + ]: 1502 : if (optcommkpt == 1) then
1687 [ + + ]: 1289 : recvcounts(:) = siz_buf * distrib%nw_mem_kptparal(:)
1688 : : else
1689 [ + + ]: 3834 : recvcounts(:) = siz_buf * distrib%nw_mem(:)
1690 : : end if
1691 : 1502 : displs(1) = 0
1692 [ + + ]: 3621 : do irank=2,nproc2
1693 : 3621 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
1694 : : end do ! irank
1695 : :
1696 [ + + + + ]: 1502 : if (optcommkpt == 1 .and. recvcounts(myproc2+1) == 0) then
1697 [ - + ]: 16 : siz_buf = siz_buf * merge(nw,distrib%nw_mem_kptparal(mod(paw_dmft%myproc,nproc_freq)+1),nproc_freq<=1)
1698 : : else
1699 : 1486 : siz_buf = recvcounts(myproc2+1)
1700 : : end if
1701 : :
1702 [ - + + - : 4506 : ABI_MALLOC(buffer,(siz_buf))
+ - + - -
+ - + ]
1703 [ - + + - : 4506 : ABI_MALLOC(buffer_tot,(recvcounts(nproc2)+displs(nproc2)))
+ - + - -
+ - + ]
1704 : :
1705 [ + + ]: 9943202 : buffer(:) = czero
1706 : :
1707 : : ibuf = 0
1708 [ + + ]: 396829 : do ifreq=1,nw
1709 [ + + ]: 395327 : if (optcommkpt == 1) then
1710 [ + + ]: 148566 : if (distrib%proct(ifreq) /= distrib%me_freq) cycle
1711 [ - + ]: 246761 : else if (optcommkpt == 0) then
1712 [ + + ]: 246761 : if (distrib%procf(ifreq) /= myproc) cycle
1713 : : end if ! optcommkpt
1714 [ + + ]: 644970 : do iatom=1,natom
1715 : 414142 : lpawu = paw_dmft%lpawu(iatom)
1716 [ + + ]: 414142 : if (lpawu == -1) cycle
1717 : 234574 : ndim = (2*lpawu+1) * nspinor
1718 [ + + ]: 1027561 : do isppol=1,nsppol
1719 [ + + ]: 2750542 : do im1=1,ndim
1720 [ + + ]: 11829240 : buffer(ibuf+1:ibuf+ndim) = oper(ifreq)%matlu(iatom)%mat(:,im1,isppol)
1721 : 2336400 : ibuf = ibuf + ndim
1722 : : end do ! im1
1723 : : end do ! isppol
1724 : : end do ! iatom
1725 : : end do ! ifreq
1726 : :
1727 [ + + ]: 1502 : if (optcommkpt == 1) comm = distrib%comm_freq
1728 [ + + ]: 1502 : if (present(master)) then
1729 [ - + ]: 106 : if (optcommkpt == 1) then
1730 : 0 : call xmpi_sum_master(buffer(:),master,distrib%comm_kpt,ierr)
1731 : : end if
1732 : 106 : call xmpi_gatherv(buffer(:),recvcounts(myproc2+1),buffer_tot(:),recvcounts(:),displs(:),master,comm,ierr)
1733 : : else
1734 [ + + ]: 1396 : if (optcommkpt == 1) then
1735 : 587 : call xmpi_sum(buffer(:),distrib%comm_kpt,ierr)
1736 : : end if
1737 : 1396 : call xmpi_allgatherv(buffer(:),recvcounts(myproc2+1),buffer_tot(:),recvcounts(:),displs(:),comm,ierr)
1738 : : end if ! present(master)
1739 : :
1740 : : ibuf = 0
1741 [ + + ]: 396829 : do ifreq=1,nw
1742 [ + + ]: 1094572 : do iatom=1,natom
1743 : 697743 : lpawu = paw_dmft%lpawu(iatom)
1744 [ + + ]: 697743 : if (lpawu == -1) cycle
1745 : 403775 : ndim = (2*lpawu+1) * nspinor
1746 [ + + ]: 1472560 : do isppol=1,nsppol
1747 [ + + ]: 4644651 : do im1=1,ndim
1748 [ + + ]: 19743540 : oper(ifreq)%matlu(iatom)%mat(:,im1,isppol) = buffer_tot(ibuf+1:ibuf+ndim)
1749 : 3946908 : ibuf = ibuf + ndim
1750 : : end do ! im1
1751 : : end do ! isppol
1752 : : end do ! iatom
1753 : : end do ! ifreq
1754 : :
1755 : : end if ! opt_ksloc
1756 : :
1757 [ - + ]: 1506 : ABI_FREE(recvcounts)
1758 [ - + ]: 1506 : ABI_FREE(displs)
1759 [ - + ]: 1506 : ABI_FREE(buffer)
1760 [ - + ]: 1506 : ABI_FREE(buffer_tot)
1761 : :
1762 : 1506 : end subroutine gather_oper
1763 : : !!***
1764 : :
1765 : : !!****f* m_oper/gather_oper_ks
1766 : : !! NAME
1767 : : !! gather_oper_ks
1768 : : !!
1769 : : !! FUNCTION
1770 : : !! For a single KS operator, performs a xmpi_sum on the frequency communicator,
1771 : : !! and a xmpi_allgatherv on the kpt-communicator
1772 : : !!
1773 : : !! INPUTS
1774 : : !! oper <type(oper_type)>= operator
1775 : : !! distrib <type(mpi_distrib_dmft_type)> = mpi related data
1776 : : !! opt_diag = 1 if the operator is diagonal in KS space, 0 (default) otherwise
1777 : : !!
1778 : : !! OUTPUT
1779 : : !!
1780 : : !! SOURCE
1781 : :
1782 : 4880 : subroutine gather_oper_ks(oper,distrib,paw_dmft,opt_diag)
1783 : :
1784 : : !Arguments ------------------------------------
1785 : : type(oper_type), intent(inout) :: oper
1786 : : type(mpi_distrib_dmft_type), intent(in) :: distrib
1787 : : type(paw_dmft_type), intent(in) :: paw_dmft
1788 : : integer, optional, intent(in) :: opt_diag
1789 : : !Local variables-------------------------------
1790 : : integer :: ib1,ibuf,ierr,ikpt,irank,isppol,mbandc
1791 : : integer :: me_kpt,nkpt,nproc,nproc_freq,nsppol,siz_buf
1792 : : logical :: diag
1793 : 4880 : integer, allocatable :: displs(:),recvcounts(:)
1794 : 4880 : complex(dp), allocatable :: buffer(:),buffer_tot(:)
1795 : : ! *********************************************************************
1796 : :
1797 : 4880 : mbandc = paw_dmft%mbandc
1798 : 4880 : me_kpt = distrib%me_kpt
1799 : 4880 : nkpt = paw_dmft%nkpt
1800 : 4880 : nproc = paw_dmft%nproc
1801 : 4880 : nsppol = paw_dmft%nsppol
1802 : :
1803 : 4880 : nproc_freq = nproc / nkpt
1804 : :
1805 : 4880 : diag = .false.
1806 [ - + ]: 4880 : if (present(opt_diag)) then
1807 [ + + ]: 4880 : if (opt_diag == 1) diag = .true.
1808 : : end if
1809 : :
1810 [ - + + - : 14640 : ABI_MALLOC(recvcounts,(nproc))
+ - + - -
+ - + ]
1811 [ - + + - : 9760 : ABI_MALLOC(displs,(nproc))
- + - + ]
1812 : :
1813 [ + + + + ]: 25542 : recvcounts(:) = distrib%nkpt_mem(:) * merge(mbandc,mbandc**2,diag)
1814 : :
1815 : 4880 : displs(1) = 0
1816 [ + + ]: 15968 : do irank=2,nproc
1817 : 15968 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
1818 : : end do ! irank
1819 : :
1820 : 4880 : siz_buf = recvcounts(me_kpt+1)
1821 [ + + ]: 4880 : if (siz_buf == 0) then
1822 : 121 : siz_buf = mbandc
1823 [ + - ]: 121 : if (nproc_freq <= 1) siz_buf = siz_buf * distrib%nkpt_mem(mod(paw_dmft%myproc,min(nkpt,nproc))+1)
1824 [ + + ]: 121 : if (.not. diag) siz_buf = siz_buf * mbandc
1825 : : end if
1826 [ - + + - : 14640 : ABI_MALLOC(buffer,(siz_buf))
+ - + - -
+ - + ]
1827 [ - + + - : 14640 : ABI_MALLOC(buffer_tot,(recvcounts(nproc)+displs(nproc)))
+ - + - -
+ - + ]
1828 : :
1829 [ + + ]: 11564 : do isppol=1,nsppol
1830 : :
1831 [ + + ]: 2975136 : buffer(:) = czero
1832 : :
1833 : : ibuf = 0
1834 [ + + ]: 55390 : do ikpt=1,nkpt
1835 [ + + ]: 48706 : if (distrib%procb(ikpt) /= me_kpt) cycle
1836 [ + + ]: 227635 : do ib1=1,mbandc
1837 [ + + ]: 248614 : if (diag) then
1838 : 8352 : ibuf = ibuf + 1
1839 : 8352 : buffer(ibuf) = oper%ks(ib1,ib1,ikpt,isppol)
1840 : : else
1841 [ + + ]: 2993264 : buffer(ibuf+1:ibuf+mbandc) = oper%ks(:,ib1,ikpt,isppol)
1842 : : ibuf = ibuf + mbandc
1843 : : end if ! diag
1844 : : end do ! ib1
1845 : : end do ! ikpt
1846 : :
1847 : 6684 : call xmpi_sum(buffer(:),distrib%comm_freq,ierr)
1848 : :
1849 : : call xmpi_allgatherv(buffer(:),recvcounts(me_kpt+1),&
1850 : 6684 : & buffer_tot(:),recvcounts(:),displs(:),distrib%comm_kpt,ierr)
1851 : :
1852 : 6684 : ibuf = 0
1853 [ + + ]: 60270 : do ikpt=1,nkpt
1854 [ + + ]: 534886 : do ib1=1,mbandc
1855 [ + + ]: 528202 : if (diag) then
1856 : 19704 : ibuf = ibuf + 1
1857 : 19704 : oper%ks(ib1,ib1,ikpt,isppol) = buffer_tot(ibuf)
1858 : : else
1859 [ + + ]: 6788816 : oper%ks(:,ib1,ikpt,isppol) = buffer_tot(ibuf+1:ibuf+mbandc)
1860 : : ibuf = ibuf + mbandc
1861 : : end if ! diag
1862 : : end do ! ib
1863 : : end do ! ikpt
1864 : :
1865 : : end do ! isppol
1866 : :
1867 [ - + ]: 4880 : ABI_FREE(buffer)
1868 [ - + ]: 4880 : ABI_FREE(buffer_tot)
1869 : 4880 : ABI_FREE(recvcounts)
1870 : 4880 : ABI_FREE(displs)
1871 : :
1872 : 4880 : end subroutine gather_oper_ks
1873 : : !!***
1874 : :
1875 [ # # # # : 0 : END MODULE m_oper
# # # # #
# # # # #
# # # # #
# ]
1876 : : !!***
|
