Branch data Line data Source code
1 : : !!****m* ABINIT/m_cplxtools
2 : : !! NAME
3 : : !! m_cplxtools
4 : : !!
5 : : !! FUNCTION
6 : : !! This module defines helper functions to operate on complex arrays (mainly used in the GW code)
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1992-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : : !!
15 : : !! NOTES
16 : : !! 1) The convention about names of interfaced routine is: cplx_<name>,
17 : : !! where <name> is equal to the name of the standard BLAS routine
18 : : !!
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_cplxtools
26 : :
27 : : use, intrinsic :: iso_c_binding
28 : : use defs_basis
29 : : use m_abicore
30 : : use m_errors
31 : : use m_xmpi
32 : :
33 : : use m_fstrings, only : toupper, sjoin
34 : :
35 : : implicit none
36 : :
37 : : private
38 : :
39 : : ! Helper functions.
40 : : public :: cplx_fromreal
41 : : public :: cplx_filter
42 : :
43 : : ! Blas1
44 : : public :: cplx_real_zdotc
45 : : public :: cplx_zaxpby
46 : :
47 : : ! Blas2
48 : : public :: cplx_zgemv
49 : :
50 : : !Blas3
51 : : public :: cplx_zgemm
52 : :
53 : : ! Helper functions for DFT calculations.
54 : : public :: cplx_box2gsph
55 : : public :: cplx_gsph2box
56 : : public :: cplx_setaug_zero
57 : : public :: cplx_setaug_zero_dpc
58 : : public :: cplx_setaug_zero_spc
59 : : public :: cplx_addtorho
60 : : public :: cplx_mat_plus_bc
61 : : !***
62 : :
63 : : ! Interfaces
64 : : interface cplx_box2gsph
65 : : module procedure cplx_box2gsph_spc
66 : : module procedure cplx_box2gsph_dpc
67 : : end interface cplx_box2gsph
68 : :
69 : : interface cplx_gsph2box
70 : : module procedure cplx_gsph2box_spc
71 : : module procedure cplx_gsph2box_dpc
72 : : end interface cplx_gsph2box
73 : :
74 : : interface cplx_setaug_zero
75 : : module procedure cplx_setaug_zero_spc
76 : : module procedure cplx_setaug_zero_dpc
77 : : end interface cplx_setaug_zero
78 : :
79 : : interface cplx_addtorho
80 : : module procedure cplx_addtorho_dpc
81 : : end interface cplx_addtorho
82 : :
83 : : interface cplx_mat_plus_bc
84 : : module procedure cplx_mat_plus_bc_spc
85 : : module procedure cplx_mat_plus_bc_dpc
86 : : end interface cplx_mat_plus_bc
87 : :
88 : : !integer,parameter,private :: MIN_SIZE = 5000
89 : :
90 : : complex(sp),private,parameter :: czero_spc = (0._sp,0._sp)
91 : : complex(sp),private,parameter :: cone_spc = (1._sp,0._sp)
92 : : !complex(sp) ,parameter :: j_spc=(0._sp,1.0_sp)
93 : :
94 : : complex(dp),private,parameter :: czero_dpc = (0._dp,0._dp)
95 : : complex(dp),private,parameter :: cone_dpc = (1._dp,0._dp)
96 : : !complex(dp) ,parameter :: j_dpc=(0._dp,1.0_dp)
97 : :
98 : : CONTAINS !========================================================================================
99 : : !!***
100 : :
101 : : !----------------------------------------------------------------------
102 : :
103 : : !!****f* m_cplxtools/cplx_fromreal
104 : : !! NAME
105 : : !! cplx_fromreal
106 : : !!
107 : : !! FUNCTION
108 : : !! Convert a real array with (real,imag) part to a complex array
109 : : !!
110 : : !! INPUTS
111 : : !! n = Specifies the number of elements in ocplx
112 : : !! ireal(2*n)=Input real array.
113 : : !!
114 : : !! OUTPUT
115 : : !! ocplx(n)=Output complex array
116 : : !!
117 : : !! SOURCE
118 : :
119 : 0 : subroutine cplx_fromreal(n,ireal,ocplx)
120 : :
121 : : !Arguments ------------------------------------
122 : : !scalars
123 : : integer,intent(in) :: n
124 : : !arrays
125 : : real(dp),intent(in) :: ireal(2,n)
126 : : complex(dp),intent(out) :: ocplx(n)
127 : :
128 : : !Local variables ------------------------------
129 : : !scalars
130 : : integer :: ii
131 : :
132 : : ! *************************************************************************
133 : :
134 : : !$OMP PARALLEL DO PRIVATE(ii)
135 [ # # ]: 0 : do ii=1,n
136 : 0 : ocplx(ii) = DCMPLX(ireal(1,ii),ireal(2,ii))
137 : : end do
138 : :
139 : 0 : end subroutine cplx_fromreal
140 : : !!***
141 : :
142 : : !----------------------------------------------------------------------
143 : :
144 : : !!****f* m_cplxtools/cplx_filter
145 : : !! NAME
146 : : !! cplx_filter
147 : : !!
148 : : !! FUNCTION
149 : : !! Set all the elements of x to zero where mask is .TRUE.
150 : : !!
151 : : !! INPUTS
152 : : !! n=Specifies the number of elements in vectors x and y.
153 : : !! mask(n)=Logical array.
154 : : !!
155 : : !! SOURCE
156 : :
157 : 0 : subroutine cplx_filter(n, x, mask)
158 : :
159 : : !Arguments ------------------------------------
160 : : !scalars
161 : : integer,intent(in) :: n
162 : : !arrays
163 : : complex(dp),intent(inout) :: x(n)
164 : : logical,intent(in) :: mask(n)
165 : :
166 : : ! *************************************************************************
167 : :
168 [ # # # # ]: 0 : where (mask)
169 : : x = czero
170 : : end where
171 : :
172 : 0 : end subroutine cplx_filter
173 : : !!***
174 : :
175 : : !----------------------------------------------------------------------
176 : :
177 : : !!****f* m_cplxtools/cplx_real_zdotc
178 : : !! NAME
179 : : !! cplx_real_zdotc
180 : : !!
181 : : !! FUNCTION
182 : : !! Perform a vector-vector operation defined as res = REAL (\Sigma (conjg(x)*y)) where x and y are n-element vectors.
183 : : !!
184 : : !! INPUTS
185 : : !! n = Specifies the number of elements in vector x and y
186 : : !! x,y = Input arrays.
187 : : !!
188 : : !! OUTPUT
189 : : !! res=Real part of the scalar product.
190 : : !!
191 : : !! SOURCE
192 : :
193 : 0 : function cplx_real_zdotc(n,x,y) result(res)
194 : :
195 : : !Arguments ------------------------------------
196 : : !scalars
197 : : integer,intent(in) :: n
198 : : !arrays
199 : : complex(dp),intent(in) :: x(n)
200 : : complex(dp),intent(in) :: y(n)
201 : : real(dp) :: res
202 : :
203 : : !Local variables-------------------------------
204 : : real(dp),external :: ddot
205 : :
206 : : ! *************************************************************************
207 : :
208 : 0 : res = ddot(2*n,x,1,y,1)
209 : :
210 : 0 : end function cplx_real_zdotc
211 : : !!***
212 : :
213 : : !----------------------------------------------------------------------
214 : :
215 : : !!****f* m_cplxtools/cplx_zaxpby
216 : : !! NAME
217 : : !! cplx_zaxpby
218 : : !!
219 : : !! FUNCTION
220 : : !! Scales two vectors, adds them to one another and stores result in the vector.
221 : : !! y := a*x + b*y
222 : : !!
223 : : !! INPUTS
224 : : !! n = the number of elements in vectors x and y.
225 : : !! a = Specifies the scalar a.
226 : : !! x = Array.
227 : : !! b = Specifies the scalar b.
228 : : !! y = Array
229 : : !!
230 : : !! OUTPUT
231 : : !! y Contains the updated vector y.
232 : : !!
233 : : !! SOURCE
234 : :
235 : 0 : subroutine cplx_zaxpby(n,a,x,b,y)
236 : :
237 : : !Arguments ------------------------------------
238 : : !scalars
239 : : integer,intent(in) :: n
240 : : complex(dp),intent(in) :: a,b
241 : : !arrays
242 : : complex(dp),intent(in) :: x(n)
243 : : complex(dp),intent(inout) :: y(n)
244 : : ! *************************************************************************
245 : :
246 : : #ifdef HAVE_LINALG_AXPBY
247 : 0 : call zaxpby(n, a, x, 1, b, y, 1)
248 : : #else
249 : : call zscal(n, b, y, 1)
250 : : call zaxpy(n, a, x, 1, y,1)
251 : : #endif
252 : :
253 : 0 : end subroutine cplx_zaxpby
254 : : !!***
255 : :
256 : : !----------------------------------------------------------------------
257 : :
258 : : !!****f* m_cplxtools/cplx_zgemv
259 : : !! NAME
260 : : !! cplx_zgemv
261 : : !!
262 : : !! FUNCTION
263 : : !! The ?gemv routines perform a matrix-vector operation defined as
264 : : !!
265 : : !! y := alpha*A*x + beta*y,
266 : : !! or
267 : : !! y := alpha*A'*x + beta*y,
268 : : !! or
269 : : !! y := alpha*conjg(A')*x + beta*y,
270 : : !!
271 : : !! where: alpha and beta are scalars, x and y are vectors, A is an m-by-n matrix.
272 : : !!
273 : : !! INPUTS
274 : : !!
275 : : !! OUTPUT
276 : : !!
277 : : !! SOURCE
278 : :
279 : 0 : subroutine cplx_zgemv(trans,nrows,ncols,mat,vec,matvec,alpha,beta)
280 : :
281 : : !Arguments ------------------------------------
282 : : !scalars
283 : : integer,intent(in) :: nrows,ncols
284 : : complex(dp),optional,intent(in) :: alpha,beta
285 : : character(len=1),intent(in) :: trans
286 : : !arrays
287 : : complex(dp),intent(in) :: mat(nrows*ncols)
288 : : complex(dp),intent(in) :: vec(*)
289 : : complex(dp),intent(inout) :: matvec(*)
290 : :
291 : : !Local variables-------------------------------
292 : : !scalars
293 : : integer :: mm,nn,kk,lda,ldb,ldc
294 : : complex(dp) :: my_alpha,my_beta
295 : : ! *************************************************************************
296 : :
297 : 0 : lda = nrows
298 : 0 : mm = nrows
299 : 0 : nn = 1
300 : 0 : kk = ncols
301 : :
302 [ # # ]: 0 : if (toupper(trans) /= 'N') then
303 : 0 : mm = ncols
304 : 0 : kk = nrows
305 : : end if
306 : :
307 : 0 : ldb = kk
308 : 0 : ldc = mm
309 : :
310 [ # # ]: 0 : my_alpha = cone_dpc; if (PRESENT(alpha)) my_alpha = alpha
311 [ # # ]: 0 : my_beta = czero_dpc; if (PRESENT(beta)) my_beta = beta
312 : :
313 : 0 : call ZGEMM(trans,"N",mm,nn,kk,my_alpha,mat,lda,vec,ldb,my_beta,matvec,ldc)
314 : :
315 : : !call ZGEMV(trans,mm,nn,my_alpha,mat,lda,vec,1,my_beta,matvec,1)
316 : :
317 : 0 : end subroutine cplx_zgemv
318 : : !!***
319 : :
320 : : !----------------------------------------------------------------------
321 : :
322 : : !!****f* m_cplxtools/cplx_zgemm
323 : : !! NAME
324 : : !! cplx_zgemm
325 : : !!
326 : : !! FUNCTION
327 : : !! The ?gemm routines perform a matrix-matrix operation with general matrices.
328 : : !! The operation is defined as C := alpha*op(A)*op(B) + beta*C,
329 : : !! where:
330 : : !!
331 : : !! op(x) is one of op(x) = x, or op(x) = x', or op(x) = conjg(x'),
332 : : !!
333 : : !! alpha and beta are scalars,
334 : : !! A, B and C are matrices:
335 : : !! op(A) is an m-by-k matrix,
336 : : !! op(B) is a k-by-n matrix,
337 : : !! C is an m-by-n matrix.
338 : : !!
339 : : !! INPUTS
340 : : !!
341 : : !! OUTPUT
342 : : !!
343 : : !! SOURCE
344 : :
345 : 0 : subroutine cplx_zgemm(transa,transb,npws,ncola,ncolb,amat,bmat,cmat,alpha,beta)
346 : :
347 : : !Arguments ------------------------------------
348 : : !scalars
349 : : integer,intent(in) :: npws,ncola,ncolb
350 : : complex(dp),optional,intent(in) :: alpha,beta
351 : : character(len=1),intent(in) :: transa,transb
352 : : !arrays
353 : : complex(dp),intent(in) :: amat(npws*ncola)
354 : : complex(dp),intent(in) :: bmat(npws*ncolb)
355 : : complex(dp),intent(inout) :: cmat(*)
356 : :
357 : : !Local variables-------------------------------
358 : : !scalars
359 : : integer :: mm,nn,kk,lda,ldb,ldc
360 : : complex(dp) :: my_alpha,my_beta
361 : : ! *************************************************************************
362 : :
363 : 0 : lda = npws
364 : 0 : ldb = npws
365 : :
366 : 0 : mm = npws
367 : 0 : nn = ncolb
368 : 0 : kk = ncola
369 : :
370 [ # # ]: 0 : if (toupper(transa) /= 'N') then
371 : 0 : mm = ncola
372 : 0 : kk = npws
373 : : end if
374 [ # # ]: 0 : if (toupper(transb) /= 'N') nn = npws
375 : :
376 : 0 : ldc = mm
377 : :
378 [ # # ]: 0 : my_alpha = cone_dpc; if (PRESENT(alpha)) my_alpha = alpha
379 [ # # ]: 0 : my_beta = czero_dpc; if (PRESENT(beta)) my_beta = beta
380 : :
381 : 0 : call ZGEMM(transa,transb,mm,nn,kk,my_alpha,amat,lda,bmat,ldb,my_beta,cmat,ldc)
382 : :
383 : 0 : end subroutine cplx_zgemm
384 : : !!***
385 : :
386 : : !----------------------------------------------------------------------
387 : :
388 : : !!****f* m_cplxtools/cplx_box2gsph_spc
389 : : !! NAME
390 : : !! cplx_box2gsph_spc
391 : : !!
392 : : !! FUNCTION
393 : : !! Transfer data from the FFT box to the G-sphere. Target SPC complex array.
394 : : !!
395 : : !! INPUTS
396 : : !! nx,ny,nz=physical dimension of the FFT box.
397 : : !! ldx,ldy,ldz=Logical dimensions of the arrays.
398 : : !! ndat=number of data in iarrbox
399 : : !! npw_k=Number of planewaves in the G-sphere.
400 : : !! kg_k(3,npw_k)=Reduced coordinates of the G-vectoes.
401 : : !! iarrbox(ldx*ldy*ldz*ndat)=Complex Input arrays on the FFT box.
402 : : !! [rscal] = Scaling factor
403 : : !!
404 : : !! OUTPUT
405 : : !! oarrsph(npw_k*ndat)=Complex Data defined on the G-sphere.
406 : : !!
407 : : !! SOURCE
408 : :
409 : 16 : subroutine cplx_box2gsph_spc(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal)
410 : :
411 : : !Arguments ------------------------------------
412 : : !scalars
413 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat
414 : : real(sp),optional,intent(in) :: rscal
415 : : !arrays
416 : : integer,intent(in) :: kg_k(3,npw_k)
417 : : complex(sp),intent(in) :: iarrbox(ldx*ldy*ldz*ndat)
418 : : complex(sp),intent(out) :: oarrsph(npw_k*ndat)
419 : :
420 : : !Local variables-------------------------------
421 : : integer :: ig,ix,iy,iz,dat,pad_sph,pad_box,ifft,ldxyz
422 : : ! *************************************************************************
423 : :
424 : 16 : ldxyz = ldx*ldy*ldz
425 [ + - ]: 16 : if (.not. PRESENT(rscal)) then
426 : : !
427 [ + + ]: 16 : if (ndat==1) then
428 : : !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
429 [ + + ]: 251248 : do ig=1,npw_k
430 [ + + ]: 251240 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
431 [ - + ]: 251240 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
432 [ + + ]: 251240 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
433 : 251240 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
434 : 251248 : oarrsph(ig) = iarrbox(ifft)
435 : : end do
436 : : else
437 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
438 [ + + ]: 40 : do dat=1,ndat
439 : 32 : pad_sph = (dat-1)*npw_k
440 : 32 : pad_box = (dat-1)*ldxyz
441 [ + + ]: 1005000 : do ig=1,npw_k
442 [ + + ]: 1004960 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
443 [ - + ]: 1004960 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
444 [ + + ]: 1004960 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
445 : 1004960 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
446 : 1004992 : oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box)
447 : : end do
448 : : end do
449 : : end if
450 : : !
451 : : else
452 [ # # ]: 0 : if (ndat==1) then
453 : : !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
454 [ # # ]: 0 : do ig=1,npw_k
455 [ # # ]: 0 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
456 [ # # ]: 0 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
457 [ # # ]: 0 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
458 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
459 : 0 : oarrsph(ig) = iarrbox(ifft) * rscal
460 : : end do
461 : : else
462 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
463 [ # # ]: 0 : do dat=1,ndat
464 : 0 : pad_sph = (dat-1)*npw_k
465 : 0 : pad_box = (dat-1)*ldxyz
466 [ # # ]: 0 : do ig=1,npw_k
467 [ # # ]: 0 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
468 [ # # ]: 0 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
469 [ # # ]: 0 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
470 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
471 : 0 : oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box) * rscal
472 : : end do
473 : : end do
474 : : end if
475 : : end if
476 : :
477 : 16 : end subroutine cplx_box2gsph_spc
478 : : !!***
479 : :
480 : : !----------------------------------------------------------------------
481 : :
482 : : !----------------------------------------------------------------------
483 : :
484 : : !!****f* m_cplxtools/cplx_box2gsph_dpc
485 : : !! NAME
486 : : !! cplx_box2gsph_dpc
487 : : !!
488 : : !! FUNCTION
489 : : !! Transfer data from the FFT box to the G-sphere. Target DP complex array.
490 : : !!
491 : : !! INPUTS
492 : : !! nx,ny,nz=physical dimension of the FFT box.
493 : : !! ldx,ldy,ldz=Logical dimensions of the arrays.
494 : : !! ndat=number of data in iarrbox
495 : : !! npw_k=Number of planewaves in the G-sphere.
496 : : !! kg_k(3,npw_k)=Reduced coordinates of the G-vectoes.
497 : : !! iarrbox(ldx*ldy*ldz*ndat)=Complex Input arrays on the FFT box.
498 : : !! [rscal] = Scaling factor
499 : : !!
500 : : !! OUTPUT
501 : : !! oarrsph(npw_k*ndat)=Complex Data defined on the G-sphere.
502 : : !!
503 : : !! SOURCE
504 : :
505 : 16 : subroutine cplx_box2gsph_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal)
506 : :
507 : : !Arguments ------------------------------------
508 : : !scalars
509 : : integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat
510 : : real(dp),optional,intent(in) :: rscal
511 : : !arrays
512 : : integer,intent(in) :: kg_k(3,npw_k)
513 : : complex(dp),intent(in) :: iarrbox(ldx*ldy*ldz*ndat)
514 : : complex(dp),intent(out) :: oarrsph(npw_k*ndat)
515 : :
516 : : !Local variables-------------------------------
517 : : !scalars
518 : : integer :: ig,ix,iy,iz,dat,pad_sph,pad_box,ifft,ldxyz
519 : :
520 : : ! *************************************************************************
521 : :
522 : 16 : ldxyz = ldx*ldy*ldz
523 [ + - ]: 16 : if (.not. PRESENT(rscal)) then
524 : : !
525 [ + + ]: 16 : if (ndat==1) then
526 : : !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
527 [ + + ]: 251248 : do ig=1,npw_k
528 [ + + ]: 251240 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
529 [ - + ]: 251240 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
530 [ + + ]: 251240 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
531 : 251240 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
532 : 251248 : oarrsph(ig) = iarrbox(ifft)
533 : : end do
534 : : else
535 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
536 [ + + ]: 40 : do dat=1,ndat
537 : 32 : pad_sph = (dat-1)*npw_k
538 : 32 : pad_box = (dat-1)*ldxyz
539 [ + + ]: 1005000 : do ig=1,npw_k
540 [ + + ]: 1004960 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
541 [ - + ]: 1004960 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
542 [ + + ]: 1004960 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
543 : 1004960 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
544 : 1004992 : oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box)
545 : : end do
546 : : end do
547 : : end if
548 : : !
549 : : else
550 [ # # ]: 0 : if (ndat==1) then
551 : : !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
552 [ # # ]: 0 : do ig=1,npw_k
553 [ # # ]: 0 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
554 [ # # ]: 0 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
555 [ # # ]: 0 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
556 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
557 : 0 : oarrsph(ig) = iarrbox(ifft) * rscal
558 : : end do
559 : : else
560 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
561 [ # # ]: 0 : do dat=1,ndat
562 : 0 : pad_sph = (dat-1)*npw_k
563 : 0 : pad_box = (dat-1)*ldxyz
564 [ # # ]: 0 : do ig=1,npw_k
565 [ # # ]: 0 : ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
566 [ # # ]: 0 : iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
567 [ # # ]: 0 : iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
568 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
569 : 0 : oarrsph(ig+pad_sph) = iarrbox(ifft+pad_box) * rscal
570 : : end do
571 : : end do
572 : : end if
573 : : end if
574 : :
575 : 16 : end subroutine cplx_box2gsph_dpc
576 : : !!***
577 : :
578 : : !----------------------------------------------------------------------
579 : :
580 : : !!****f* m_cplxtools/cplx_gsph2box_spc
581 : : !! NAME
582 : : !! cplx_gsph2box_spc
583 : : !!
584 : : !! FUNCTION
585 : : !! Array iarrsph is defined in sphere with npw points. Insert iarrsph inside box
586 : : !! of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.
587 : : !! targer: SPC complex arrays
588 : : !!
589 : : !! INPUTS
590 : : !! iarrsph(2,npw*ndat)= contains values for npw G vectors in basis sphere
591 : : !! ndat=number of FFT to perform.
592 : : !! npw=number of G vectors in basis at this k point
593 : : !! oarrbox(2,ldx*ldy*ldz*ndat) = fft box
594 : : !! nx,ny,nz=physical dimension of the box (oarrbox)
595 : : !! ldx,ldy,ldz=memory dimension of oarrbox
596 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
597 : : !! istwf_k=option parameter that describes the storage of wfs
598 : : !!
599 : : !! OUTPUT
600 : : !! oarrbox(ldx*ldy*ldz*ndat)
601 : : !!
602 : : !! NOTES
603 : : !! If istwf_k differs from 1, then special storage modes must be taken
604 : : !! into account, for symmetric wavefunctions coming from k=(0 0 0) or other
605 : : !! special k points.
606 : : !!
607 : : !! SOURCE
608 : :
609 : 14 : subroutine cplx_gsph2box_spc(nx,ny,nz,ldx,ldy,ldz,ndat,npw,istwf_k,kg_k,iarrsph,oarrbox)
610 : :
611 : : !Arguments ------------------------------------
612 : : !scalars
613 : : integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw
614 : : !arrays
615 : : integer,intent(in) :: kg_k(3,npw)
616 : : complex(sp),intent(in) :: iarrsph(npw*ndat)
617 : : complex(sp),intent(out) :: oarrbox(ldx*ldy*ldz*ndat)
618 : :
619 : : !Local variables-------------------------------
620 : : !scalars
621 : : integer,parameter :: me_g0=1
622 : : integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz
623 : : !character(len=500) :: msg
624 : : !arrays
625 : 14 : integer,allocatable :: ixinver(:),iyinver(:),izinver(:)
626 : : ! *************************************************************************
627 : :
628 : : !In the case of special k-points, invariant under time-reversal,
629 : : !but not Gamma, initialize the inverse coordinates
630 : : !Remember indeed that
631 : : !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
632 : : !and therefore:
633 : : !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
634 [ + - ]: 14 : if (istwf_k>=2) then
635 [ - + + - : 42 : ABI_MALLOC(ixinver,(nx))
+ - + - -
+ - + ]
636 [ - + + - : 42 : ABI_MALLOC(iyinver,(ny))
+ - + - -
+ - + ]
637 [ - + + - : 42 : ABI_MALLOC(izinver,(nz))
+ - + - -
+ - + ]
638 [ + + ]: 14 : if ( ANY(istwf_k==(/2,4,6,8/)) ) then
639 : 6 : ixinver(1)=1
640 [ + + ]: 600 : do ix=2,nx
641 : 600 : ixinver(ix)=nx+2-ix
642 : : end do
643 : : else
644 [ + + ]: 808 : do ix=1,nx
645 : 808 : ixinver(ix)=nx+1-ix
646 : : end do
647 : : end if
648 [ + + ]: 14 : if (istwf_k>=2 .and. istwf_k<=5) then
649 : 6 : iyinver(1)=1
650 [ + + ]: 600 : do iy=2,ny
651 : 600 : iyinver(iy)=ny+2-iy
652 : : end do
653 : : else
654 [ + + ]: 808 : do iy=1,ny
655 : 808 : iyinver(iy)=ny+1-iy
656 : : end do
657 : : end if
658 [ + + + + ]: 14 : if ( ANY(istwf_k==(/2,3,6,7/)) ) then
659 : 6 : izinver(1)=1
660 [ + + ]: 600 : do iz=2,nz
661 : 600 : izinver(iz)=nz+2-iz
662 : : end do
663 : : else
664 [ + + ]: 808 : do iz=1,nz
665 : 808 : izinver(iz)=nz+1-iz
666 : : end do
667 : : end if
668 : : end if
669 : :
670 : 14 : ldxyz = ldx*ldy*ldz
671 : :
672 [ - + ]: 14 : if (istwf_k==1) then
673 : :
674 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
675 [ # # ]: 0 : do dat=1,ndat
676 : 0 : pad_sph = (dat-1)*npw
677 : 0 : pad_box = (dat-1)*ldxyz
678 [ # # ]: 0 : oarrbox(1+pad_box:ldxyz+pad_box) = czero_spc ! zero the sub-array
679 [ # # ]: 0 : do ipw=1,npw
680 [ # # ]: 0 : ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1
681 [ # # ]: 0 : iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1
682 [ # # ]: 0 : iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1
683 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
684 : : #if defined __INTEL_COMPILER && defined HAVE_OPENMP
685 : : if (ifft==0) then
686 : : ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
687 : : end if
688 : : #endif
689 : 0 : oarrbox(ifft+pad_box) = iarrsph(ipw+pad_sph)
690 : : end do
691 : : end do
692 : :
693 [ + - ]: 14 : else if (istwf_k>=2) then
694 : : !
695 : 14 : npwmin=1
696 [ - + ]: 14 : if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed
697 [ # # ]: 0 : do dat=1,ndat
698 : 0 : pad_sph = (dat-1)*npw
699 : 0 : pad_box = (dat-1)*ldxyz
700 : 0 : oarrbox(1+pad_box) = REAL(iarrsph(1+pad_sph))
701 : : end do
702 : : npwmin=2
703 : : end if
704 : :
705 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,ixinv,iy,iyinv,iz,izinv,ifft)
706 [ + + ]: 49 : do dat=1,ndat
707 : 35 : pad_sph = (dat-1)*npw
708 : 35 : pad_box = (dat-1)*ldxyz
709 [ + + ]: 35000035 : oarrbox(npwmin+pad_box:ldxyz+pad_box) = czero_spc
710 [ + + ]: 1099574 : do ipw=npwmin,npw
711 [ + + ]: 1099525 : ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1
712 [ - + ]: 1099525 : iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1
713 [ + + ]: 1099525 : iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1
714 : 1099525 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
715 : : ! Construct the coordinates of -k-G
716 : 1099525 : ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz)
717 : 1099525 : ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy
718 : : #if defined __INTEL_COMPILER && defined HAVE_OPENMP
719 : : if (ifft==0 .or. ifft_inv==0) then
720 : : ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
721 : : end if
722 : : #endif
723 : 1099525 : oarrbox(ifft +pad_box) = iarrsph(ipw+pad_sph)
724 : 1099560 : oarrbox(ifft_inv+pad_box) = CONJG(iarrsph(ipw+pad_sph))
725 : : end do
726 : : end do
727 : : !
728 : : else
729 : 0 : ABI_ERROR("Wrong istwfk")
730 : : end if
731 : :
732 [ + - ]: 14 : if (istwf_k>=2) then
733 [ - + ]: 14 : ABI_FREE(ixinver)
734 [ - + ]: 14 : ABI_FREE(iyinver)
735 [ - + ]: 14 : ABI_FREE(izinver)
736 : : end if
737 : :
738 [ - + - + : 14 : end subroutine cplx_gsph2box_spc
- + ]
739 : : !!***
740 : :
741 : : !----------------------------------------------------------------------
742 : :
743 : : !!****f* m_cplxtools/cplx_gsph2box_dpc
744 : : !! NAME
745 : : !! cplx_gsph2box_dpc
746 : : !!
747 : : !! FUNCTION
748 : : !! Array iarrsph is defined in sphere with npw points. Insert iarrsph inside box
749 : : !! of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.
750 : : !! targer: DP complex arrays
751 : : !!
752 : : !! INPUTS
753 : : !! iarrsph(2,npw*ndat)= contains values for npw G vectors in basis sphere
754 : : !! ndat=number of FFT to perform.
755 : : !! npw=number of G vectors in basis at this k point
756 : : !! oarrbox(2,ldx*ldy*ldz*ndat) = fft box
757 : : !! nx,ny,nz=physical dimension of the box (oarrbox)
758 : : !! ldx,ldy,ldz=memory dimension of oarrbox
759 : : !! kg_k(3,npw)=integer coordinates of G vectors in basis sphere
760 : : !! istwf_k=option parameter that describes the storage of wfs
761 : : !!
762 : : !! OUTPUT
763 : : !! oarrbox(ldx*ldy*ldz*ndat)
764 : : !!
765 : : !! NOTES
766 : : !! If istwf_k differs from 1, then special storage modes must be taken
767 : : !! into account, for symmetric wavefunctions coming from k=(0 0 0) or other
768 : : !! special k points.
769 : : !!
770 : : !! SOURCE
771 : :
772 : 2016 : subroutine cplx_gsph2box_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,npw,istwf_k,kg_k,iarrsph,oarrbox)
773 : :
774 : : !Arguments ------------------------------------
775 : : !scalars
776 : : integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw
777 : : !arrays
778 : : integer,intent(in) :: kg_k(3,npw)
779 : : complex(dp),intent(in) :: iarrsph(npw*ndat)
780 : : complex(dp),intent(out) :: oarrbox(ldx*ldy*ldz*ndat)
781 : :
782 : : !Local variables-------------------------------
783 : : !scalars
784 : : integer,parameter :: me_g0=1
785 : : integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz
786 : : !character(len=500) :: msg
787 : : !arrays
788 : 2016 : integer,allocatable :: ixinver(:),iyinver(:),izinver(:)
789 : : ! *************************************************************************
790 : :
791 : : !In the case of special k-points, invariant under time-reversal,
792 : : !but not Gamma, initialize the inverse coordinates
793 : : !Remember indeed that
794 : : !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
795 : : !and therefore:
796 : : !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
797 [ + - ]: 2016 : if (istwf_k>=2) then
798 [ - + + - : 6048 : ABI_MALLOC(ixinver,(nx))
+ - + - -
+ - + ]
799 [ - + + - : 6048 : ABI_MALLOC(iyinver,(ny))
+ - + - -
+ - + ]
800 [ - + + - : 6048 : ABI_MALLOC(izinver,(nz))
+ - + - -
+ - + ]
801 [ + + ]: 2016 : if ( ANY(istwf_k==(/2,4,6,8/)) ) then
802 : 330 : ixinver(1)=1
803 [ + + ]: 7944 : do ix=2,nx
804 : 7944 : ixinver(ix)=nx+2-ix
805 : : end do
806 : : else
807 [ + + ]: 40389 : do ix=1,nx
808 : 40389 : ixinver(ix)=nx+1-ix
809 : : end do
810 : : end if
811 [ + + ]: 2016 : if (istwf_k>=2 .and. istwf_k<=5) then
812 : 1338 : iyinver(1)=1
813 [ + + ]: 30705 : do iy=2,ny
814 : 30705 : iyinver(iy)=ny+2-iy
815 : : end do
816 : : else
817 [ + + ]: 16620 : do iy=1,ny
818 : 16620 : iyinver(iy)=ny+1-iy
819 : : end do
820 : : end if
821 [ + + + + ]: 2016 : if ( ANY(istwf_k==(/2,3,6,7/)) ) then
822 : 1208 : izinver(1)=1
823 [ + + ]: 28626 : do iz=2,nz
824 : 28626 : izinver(iz)=nz+2-iz
825 : : end do
826 : : else
827 [ + + ]: 21576 : do iz=1,nz
828 : 21576 : izinver(iz)=nz+1-iz
829 : : end do
830 : : end if
831 : : end if
832 : :
833 : 2016 : ldxyz = ldx*ldy*ldz
834 : :
835 [ - + ]: 2016 : if (istwf_k==1) then
836 : :
837 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
838 [ # # ]: 0 : do dat=1,ndat
839 : 0 : pad_sph = (dat-1)*npw
840 : 0 : pad_box = (dat-1)*ldxyz
841 [ # # ]: 0 : oarrbox(1+pad_box:ldxyz+pad_box) = czero_dpc ! zero the sub-array
842 [ # # ]: 0 : do ipw=1,npw
843 [ # # ]: 0 : ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1
844 [ # # ]: 0 : iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1
845 [ # # ]: 0 : iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1
846 : 0 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
847 : : #if defined __INTEL_COMPILER && defined HAVE_OPENMP
848 : : if (ifft==0) then
849 : : ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
850 : : end if
851 : : #endif
852 : 0 : oarrbox(ifft+pad_box) = iarrsph(ipw+pad_sph)
853 : : end do
854 : : end do
855 : :
856 [ + - ]: 2016 : else if (istwf_k>=2) then
857 : : !
858 : 2016 : npwmin=1
859 [ - + ]: 2016 : if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed
860 [ # # ]: 0 : do dat=1,ndat
861 : 0 : pad_sph = (dat-1)*npw
862 : 0 : pad_box = (dat-1)*ldxyz
863 : 0 : oarrbox(1+pad_box) = REAL(iarrsph(1+pad_sph))
864 : : end do
865 : : npwmin=2
866 : : end if
867 : :
868 : : !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,ixinv,iy,iyinv,iz,izinv,ifft)
869 [ + + ]: 4053 : do dat=1,ndat
870 : 2037 : pad_sph = (dat-1)*npw
871 : 2037 : pad_box = (dat-1)*ldxyz
872 [ + + ]: 62496801 : oarrbox(npwmin+pad_box:ldxyz+pad_box) = czero_dpc
873 [ + + ]: 2150498 : do ipw=npwmin,npw
874 [ + + ]: 2146445 : ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1
875 [ - + ]: 2146445 : iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1
876 [ + + ]: 2146445 : iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1
877 : 2146445 : ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
878 : : ! Construct the coordinates of -k-G
879 : 2146445 : ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz)
880 : 2146445 : ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy
881 : : #if defined __INTEL_COMPILER && defined HAVE_OPENMP
882 : : if (ifft==0 .or. ifft_inv==0) then
883 : : ABI_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
884 : : end if
885 : : #endif
886 : 2146445 : oarrbox(ifft +pad_box) = iarrsph(ipw+pad_sph)
887 : 2148482 : oarrbox(ifft_inv+pad_box) = DCONJG(iarrsph(ipw+pad_sph))
888 : : end do
889 : : end do
890 : : !
891 : : else
892 : 0 : ABI_ERROR("Wrong istwfk")
893 : : end if
894 : :
895 [ + - ]: 2016 : if (istwf_k>=2) then
896 [ - + ]: 2016 : ABI_FREE(ixinver)
897 [ - + ]: 2016 : ABI_FREE(iyinver)
898 [ - + ]: 2016 : ABI_FREE(izinver)
899 : : end if
900 : :
901 [ - + - + : 2016 : end subroutine cplx_gsph2box_dpc
- + ]
902 : : !!***
903 : :
904 : : !----------------------------------------------------------------------
905 : :
906 : : !!****f* m_cplxtools/cplx_setaug_zero_spc
907 : : !! NAME
908 : : !! cplx_setaug_zero_spc
909 : : !!
910 : : !! FUNCTION
911 : : !! Set to zero all elements of the array that are not in the FFT box.
912 : : !!
913 : : !! INPUTS
914 : : !! nx,ny,nz=physical dimensions of the FFT box
915 : : !! ldx,ldy,ldx=memory dimension of arr
916 : : !! ndat=number of FFTs
917 : : !!
918 : : !! SIDE EFFECT
919 : : !! arr(ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero
920 : : !!
921 : : !! SOURCE
922 : :
923 : 24 : subroutine cplx_setaug_zero_spc(nx,ny,nz,ldx,ldy,ldz,ndat,arr)
924 : :
925 : : !Arguments ------------------------------------
926 : : !scalars
927 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
928 : : !arrays
929 : : complex(sp),intent(inout) :: arr(ldx,ldy,ldz*ndat)
930 : :
931 : : !Local variables-------------------------------
932 : : integer :: iy,iz,dat,padat
933 : : ! *************************************************************************
934 : :
935 [ + + ]: 24 : if (nx /= ldx) then
936 [ + + ]: 656 : do iz=1,ldz*ndat
937 [ + + ]: 13026 : do iy=1,ldy
938 [ + + ]: 32940 : arr(nx+1:ldx,iy,iz) = czero_spc
939 : : end do
940 : : end do
941 : : end if
942 : :
943 [ + + ]: 24 : if (ny /= ldy) then
944 [ + + ]: 502 : do iz=1,ldz*ndat
945 [ + + + + ]: 13482 : arr(:,ny+1:ldy,iz) = czero_spc
946 : : end do
947 : : end if
948 : :
949 [ + + ]: 24 : if (nz /= ldz) then
950 [ + + ]: 42 : do dat=1,ndat
951 : 30 : padat = ldz*(dat-1)
952 [ + + ]: 92 : do iz=nz+1,ldz
953 [ + + + + ]: 15160 : arr(:,:,iz+padat) = czero_spc
954 : : end do
955 : : end do
956 : : end if
957 : :
958 : 24 : end subroutine cplx_setaug_zero_spc
959 : : !!***
960 : :
961 : : !----------------------------------------------------------------------
962 : :
963 : : !!****f* m_cplxtools/cplx_setaug_zero_dpc
964 : : !! NAME
965 : : !! cplx_setaug_zero_dpc
966 : : !!
967 : : !! FUNCTION
968 : : !! Set to zero all elements of the array that are not in the FFT box.
969 : : !!
970 : : !! INPUTS
971 : : !! nx,ny,nz=physical dimensions of the FFT box
972 : : !! ldx,ldy,ldx=memory dimension of arr
973 : : !! ndat=number of FFTs
974 : : !!
975 : : !! SIDE EFFECT
976 : : !! arr(ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero
977 : : !!
978 : : !! SOURCE
979 : :
980 : 24 : subroutine cplx_setaug_zero_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,arr)
981 : :
982 : : !Arguments ------------------------------------
983 : : !scalars
984 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
985 : : !arrays
986 : : complex(dp),intent(inout) :: arr(ldx,ldy,ldz*ndat)
987 : :
988 : : !Local variables-------------------------------
989 : : integer :: iy,iz,dat,padat
990 : : ! *************************************************************************
991 : :
992 [ + + ]: 24 : if (nx /= ldx) then
993 [ + + ]: 656 : do iz=1,ldz*ndat
994 [ + + ]: 13026 : do iy=1,ldy
995 [ + + ]: 32940 : arr(nx+1:ldx,iy,iz) = czero_dpc
996 : : end do
997 : : end do
998 : : end if
999 : :
1000 [ + + ]: 24 : if (ny /= ldy) then
1001 [ + + ]: 502 : do iz=1,ldz*ndat
1002 [ + + + + ]: 13482 : arr(:,ny+1:ldy,iz) = czero_dpc
1003 : : end do
1004 : : end if
1005 : :
1006 [ + + ]: 24 : if (nz /= ldz) then
1007 [ + + ]: 42 : do dat=1,ndat
1008 : 30 : padat = ldz*(dat-1)
1009 [ + + ]: 92 : do iz=nz+1,ldz
1010 [ + + + + ]: 15160 : arr(:,:,iz+padat) = czero_dpc
1011 : : end do
1012 : : end do
1013 : : end if
1014 : :
1015 : 24 : end subroutine cplx_setaug_zero_dpc
1016 : : !!***
1017 : :
1018 : : !----------------------------------------------------------------------
1019 : :
1020 : : !!****f* m_cplxtools/cplx_addtorho_dpc
1021 : : !! NAME
1022 : : !! cplx_addtorho_dpc
1023 : : !!
1024 : : !! FUNCTION
1025 : : !! Add |ur|**2 to the ground-states density rho.
1026 : : !! rho = rho + weight_r * |ur|**2
1027 : : !!
1028 : : !! INPUTS
1029 : : !! nx,ny,nz=physical dimension of the FFT box.
1030 : : !! ldx,ldy,ldz=leading dimensions of the arrays.
1031 : : !! ndat=number of contributions to accumulate.
1032 : : !! weight_r=weight used for the accumulation of the density in real space
1033 : : !! ur(ldx,ldy,ldz*ndat)=wavefunctions in real space
1034 : : !!
1035 : : !! SIDE EFFECTS
1036 : : !! rho(ldx,ldy,ldz) = contains the input density at input,
1037 : : !! modified in input with the contribution gived by ur.
1038 : : !!
1039 : : !! SOURCE
1040 : :
1041 : 0 : subroutine cplx_addtorho_dpc(nx,ny,nz,ldx,ldy,ldz,ndat,weight_r,ur,rho)
1042 : :
1043 : : !Arguments ------------------------------------
1044 : : !scalars
1045 : : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
1046 : : real(dp),intent(in) :: weight_r
1047 : : !arrays
1048 : : complex(dp),intent(in) :: ur(ldx*ldy*ldz*ndat)
1049 : : real(dp),intent(inout) :: rho(ldx*ldy*ldz)
1050 : :
1051 : : !Local variables-------------------------------
1052 : : !scalars
1053 : : integer :: ix,iy,iz,dat,ifft,ldxyz,pad_box,padz,pady
1054 : : ! *************************************************************************
1055 : :
1056 : 0 : ldxyz = ldx*ldy*ldz
1057 : :
1058 [ # # ]: 0 : if (ndat==1) then
1059 : : !$OMP PARALLEL DO PRIVATE(padz, pady, ifft)
1060 [ # # ]: 0 : do iz=1,nz
1061 : 0 : padz = (iz-1)*ldx*ldy
1062 [ # # ]: 0 : do iy=1,ny
1063 : 0 : pady = (iy-1)*ldx
1064 [ # # ]: 0 : do ix=1,nx
1065 : 0 : ifft = ix + pady + padz
1066 : 0 : rho(ifft) = rho(ifft) + weight_r * (REAL(ur(ifft))**2 + AIMAG(ur(ifft))**2)
1067 : : end do
1068 : : end do
1069 : : end do
1070 : :
1071 : : else
1072 : : ! It would be nice to use $OMP PARALLEL DO REDUCTION(+:rho)
1073 : : ! but it's risky as the private rho is allocated on the stack of the thread.
1074 : : !$OMP PARALLEL PRIVATE(pad_box, padz, pady, ifft)
1075 [ # # ]: 0 : do dat=1,ndat
1076 : 0 : pad_box = (dat-1)*ldxyz
1077 : : !$OMP DO
1078 [ # # ]: 0 : do iz=1,nz
1079 : 0 : padz = (iz-1)*ldx*ldy
1080 [ # # ]: 0 : do iy=1,ny
1081 : 0 : pady = (iy-1)*ldx
1082 [ # # ]: 0 : do ix=1,nx
1083 : 0 : ifft = ix + pady + padz
1084 : 0 : rho(ifft) = rho(ifft) + weight_r * (REAL(ur(ifft+pad_box)**2 + AIMAG(ur(ifft+pad_box))**2))
1085 : : end do
1086 : : end do
1087 : : end do
1088 : : !$OMP END DO NOWAIT
1089 : : end do
1090 : : !$OMP END PARALLEL
1091 : : end if
1092 : :
1093 : 0 : end subroutine cplx_addtorho_dpc
1094 : : !!***
1095 : :
1096 : : !----------------------------------------------------------------------
1097 : :
1098 : : !!****f* m_cplxtools/cplx_mat_plus_bc_spc
1099 : : !! NAME
1100 : : !! cplx_mat_plus_bc_spc
1101 : : !!
1102 : : !! FUNCTION
1103 : : !! Compute a = a + fact * b * c
1104 : : !!
1105 : : !! SOURCE
1106 : :
1107 : 0 : subroutine cplx_mat_plus_bc_spc(bufsize, aa, real_fact, bmode, bb, cc, gpu_option)
1108 : :
1109 : : !Arguments ------------------------------------
1110 : : !scalars
1111 : : integer(c_size_t),intent(in) :: bufsize
1112 : : integer,intent(in) :: gpu_option
1113 : : real(sp),intent(in) :: real_fact
1114 : : character(len=1),intent(in) :: bmode
1115 : : !arrays
1116 : : complex(sp),intent(inout) :: aa(bufsize)
1117 : : complex(sp),intent(in) :: bb(bufsize), cc(bufsize)
1118 : :
1119 : : !Local variables-------------------------------
1120 : : integer(c_size_t) :: ii
1121 : : ! *************************************************************************
1122 : :
1123 [ # # ]: 0 : if (gpu_option == ABI_GPU_DISABLED) then
1124 : : select case (bmode)
1125 : : case ("N")
1126 : : !$OMP PARALLEL DO
1127 [ # # ]: 0 : do ii=1,bufsize
1128 : 0 : aa(ii) = aa(ii) + real_fact * bb(ii) * cc(ii)
1129 : : end do
1130 : : case ("C")
1131 : : !$OMP PARALLEL DO
1132 [ # # ]: 0 : do ii=1,bufsize
1133 : 0 : aa(ii) = aa(ii) + real_fact * conjg(bb(ii)) * cc(ii)
1134 : : end do
1135 : : case default
1136 [ # # # ]: 0 : ABI_ERROR(sjoin("Invalid bmode:", bmode))
1137 : : end select
1138 : :
1139 : : else
1140 : : #ifdef HAVE_OPENMP_OFFLOAD
1141 : : select case (bmode)
1142 : : case ("N")
1143 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO MAP(to:aa, bb, cc)
1144 : : do ii=1,bufsize
1145 : : aa(ii) = aa(ii) + real_fact * bb(ii) * cc(ii)
1146 : : end do
1147 : : case ("C")
1148 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO MAP(to:aa, bb, cc)
1149 : : do ii=1,bufsize
1150 : : aa(ii) = aa(ii) + real_fact * conjg(bb(ii)) * cc(ii)
1151 : : end do
1152 : : case default
1153 : : ABI_ERROR(sjoin("Invalid bmode:", bmode))
1154 : : end select
1155 : : #endif
1156 : : end if
1157 : :
1158 : 0 : end subroutine cplx_mat_plus_bc_spc
1159 : : !!***
1160 : :
1161 : : !!****f* m_cplxtools/cplx_mat_plus_bc_dpc
1162 : : !! NAME
1163 : : !! cplx_mat_plus_bc_dpc
1164 : : !!
1165 : : !! FUNCTION
1166 : : !! Compute a = a + fact * b * c
1167 : : !!
1168 : : !! SOURCE
1169 : :
1170 : 0 : subroutine cplx_mat_plus_bc_dpc(bufsize, aa, real_fact, bmode, bb, cc, gpu_option)
1171 : :
1172 : : !Arguments ------------------------------------
1173 : : !scalars
1174 : : integer(c_size_t),intent(in) :: bufsize
1175 : : integer,intent(in) :: gpu_option
1176 : : real(dp),intent(in) :: real_fact
1177 : : character(len=1),intent(in) :: bmode
1178 : : !arrays
1179 : : complex(dp),intent(inout) :: aa(bufsize)
1180 : : complex(dp),intent(in) :: bb(bufsize), cc(bufsize)
1181 : :
1182 : : !Local variables-------------------------------
1183 : : integer(c_size_t) :: ii
1184 : : ! *************************************************************************
1185 : :
1186 [ # # ]: 0 : if (gpu_option == ABI_GPU_DISABLED) then
1187 : : ! CPU version
1188 : : select case (bmode)
1189 : : case ("N")
1190 : : !$OMP PARALLEL DO
1191 [ # # ]: 0 : do ii=1,bufsize
1192 : 0 : aa(ii) = aa(ii) + real_fact * bb(ii) * cc(ii)
1193 : : end do
1194 : : case ("C")
1195 : : !$OMP PARALLEL DO
1196 [ # # ]: 0 : do ii=1,bufsize
1197 : 0 : aa(ii) = aa(ii) + real_fact * conjg(bb(ii)) * cc(ii)
1198 : : end do
1199 : : case default
1200 [ # # # ]: 0 : ABI_ERROR(sjoin("Invalid bmode:", bmode))
1201 : : end select
1202 : :
1203 : : else
1204 : : ! GPU version
1205 : : #ifdef HAVE_OPENMP_OFFLOAD
1206 : : select case (bmode)
1207 : : case ("N")
1208 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO MAP(to:aa, bb, cc)
1209 : : do ii=1,bufsize
1210 : : aa(ii) = aa(ii) + real_fact * bb(ii) * cc(ii)
1211 : : end do
1212 : : case ("C")
1213 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO MAP(to:aa, bb, cc)
1214 : : do ii=1,bufsize
1215 : : aa(ii) = aa(ii) + real_fact * conjg(bb(ii)) * cc(ii)
1216 : : end do
1217 : : case default
1218 : : ABI_ERROR(sjoin("Invalid bmode:", bmode))
1219 : : end select
1220 : : #endif
1221 : : end if
1222 : :
1223 : 0 : end subroutine cplx_mat_plus_bc_dpc
1224 : : !!***
1225 : :
1226 : : end module m_cplxtools
1227 : : !!***
|
