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