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