Line data Source code
1 : !{\src2tex{textfont=tt}}
2 : !!****f* m_abi_linalg/abi_xhegv
3 : !! NAME
4 : !! abi_xhegv
5 : !!
6 : !! FUNCTION
7 : !! abi_xhegv is the generic function that compute
8 : !! all eigenvalues and, optionally, eigenvectors of a
9 : !! generalized symmetric-definite eigenproblem, of the form
10 : !! A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
11 : !! Here A and B are assumed to be symmetric (or hermitian) and
12 : !! B is also positive definite.
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2001-2026 ABINIT group (LNguyen,FDahm,MT)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !!
20 : !! SOURCE
21 :
22 : !!***
23 :
24 : !!****f* m_abi_linalg/abi_dhegv
25 : !! NAME
26 : !! abi_dhegv
27 : !!
28 : !! FUNCTION
29 : !!
30 : !! INPUTS
31 : !!
32 : !! SOURCE
33 : !!
34 :
35 28915 : subroutine abi_dhegv(itype,jobz,uplo,n,a,lda,b,ldb,w, &
36 : & x_cplx,istwf_k,timopt,tim_xeigen,&
37 : & use_slk,use_gpu_elpa,use_gpu_magma)
38 :
39 : !Arguments ------------------------------------
40 : integer, intent(in) :: itype
41 : character(len=1), intent(in) :: jobz
42 : character(len=1), intent(in) :: uplo
43 : integer, intent(in) :: n,lda,ldb
44 : real(dp), intent(inout) :: a(n,*) ! FIXME should be cplex * lda
45 : real(dp), intent(inout) :: b(n,*)
46 : real(dp),target,intent(out) :: w(n)
47 : integer, optional, intent(in) :: x_cplx
48 : integer, optional, intent(in) :: istwf_k
49 : integer, optional, intent(in) :: timopt,tim_xeigen
50 : integer, optional, intent(in) :: use_gpu_elpa,use_gpu_magma,use_slk
51 :
52 : !Local variables ------------------------------------
53 : integer :: cplx_,istwf_k_,use_gpu_elpa_,use_gpu_magma_,use_slk_
54 : integer :: info
55 : real(dp) :: tsec(2)
56 :
57 : ! *********************************************************************
58 :
59 28915 : ABI_CHECK(lapack_full_storage,"BUG(1) in abi_dhegv (storage)!")
60 28915 : ABI_CHECK(lapack_double_precision,"BUG(2) in abi_dhegv (precision)!")
61 28915 : ABI_CHECK(n<=eigen_d_maxsize,"BUG(3) in abi_dhegv (maxsize)!")
62 :
63 28915 : if (present(tim_xeigen).and.present(timopt)) then
64 28915 : if(abs(timopt)==3) call timab(tim_xeigen,1,tsec)
65 : end if
66 :
67 28915 : cplx_=1 ; if(present(x_cplx)) cplx_ = x_cplx
68 28915 : use_gpu_magma_=0; if (present(use_gpu_magma)) use_gpu_magma_=use_gpu_magma
69 28915 : use_slk_=0; if (present(use_slk)) use_slk_=use_slk
70 28915 : istwf_k_=1; if (present(istwf_k)) istwf_k_=istwf_k
71 28915 : use_gpu_elpa_=0
72 : #ifdef HAVE_LINALG_ELPA
73 : if (present(use_gpu_elpa)) use_gpu_elpa_=use_gpu_elpa
74 : #endif
75 :
76 :
77 : !===== MAGMA
78 28915 : if (ABI_LINALG_MAGMA_ISON.and.use_gpu_magma_==1) then
79 : #if defined HAVE_LINALG_MAGMA
80 : ABI_CHECK((lapack_divide_conquer),"BUG(4) in abi_dhegv (d&c)!")
81 : if (cplx_ == 2) then
82 : call magmaf_zhegvd(itype,jobz,uplo,n,a,lda,b,ldb,w,eigen_z_work,eigen_z_lwork, &
83 : & eigen_z_rwork,eigen_z_lrwork,eigen_iwork,eigen_liwork,info)
84 : else
85 : call magmaf_dsygvd(jobz,uplo,n,a,lda,w,eigen_d_work,eigen_d_lwork,&
86 : & eigen_iwork,eigen_liwork,info)
87 : endif
88 : #endif
89 :
90 : !===== SCALAPACK
91 28915 : else if (ABI_LINALG_SCALAPACK_ISON.and.use_slk_==1.and.n>slk_minsize) then
92 : #if defined HAVE_LINALG_SCALAPACK
93 : ABI_CHECK(present(x_cplx),"BUG(5) in abi_dhegv (x_cplx)!")
94 : call compute_eigen2(slk_communicator,slk_processor,cplx_,n,n,a,b,w,istwf_k_,&
95 : & use_gpu_elpa=use_gpu_elpa_)
96 : info = 0 ! This is to avoid unwanted warning but it's not clean
97 : #endif
98 :
99 : !===== PLASMA
100 : !FDahm & LNGuyen (November 2012) :
101 : ! In Plasma v 2.4.6, eigen routines support only
102 : ! the eigenvalues computation (jobz=N) and not the
103 : ! full eigenvectors bases determination (jobz=V)
104 28915 : else if (ABI_LINALG_PLASMA_ISON.and.LSAME(jobz,'N')) then
105 : #if defined HAVE_LINALG_PLASMA
106 : if (cplx_ == 2) then
107 : call PLASMA_Alloc_Workspace_zhegv(n,n,plasma_work,info)
108 : info = PLASMA_zhegv_c(itype,jobz_plasma(jobz),uplo_plasma(uplo),n,c_loc(a),lda,&
109 : & c_loc(b),ldb,c_loc(w),plasma_work,c_loc(eigen_z_work),n)
110 : else
111 : call PLASMA_Alloc_Workspace_dsygv(n,n,plasma_work,info)
112 : info = PLASMA_dsygv_c(itype,jobz_plasma(jobz),uplo_plasma(uplo),n,c_loc(a),lda,&
113 : & c_loc(b),ldb,c_loc(w),plasma_work,c_loc(eigen_d_work),n)
114 : endif
115 : call PLASMA_Dealloc_handle(plasma_work,info)
116 : #endif
117 :
118 : !===== LAPACK
119 : else
120 28915 : if (cplx_ == 2) then
121 27103 : call zhegv(itype,jobz,uplo,n,a,lda,b,ldb,w,eigen_z_work,eigen_z_lwork,eigen_z_rwork,info)
122 : else
123 1812 : call dsygv(itype,jobz,uplo,n,a,lda,b,ldb,w,eigen_d_work,eigen_d_lwork,info)
124 : endif
125 : end if
126 :
127 28915 : if (present(tim_xeigen).and.present(timopt)) then
128 28915 : if(abs(timopt)==3) call timab(tim_xeigen,2,tsec)
129 : end if
130 :
131 28915 : ABI_CHECK(info==0,"abi_dhegv returned info!=0!")
132 :
133 : #ifndef HAVE_LINALG_ELPA
134 : ABI_UNUSED(use_gpu_elpa)
135 : #endif
136 :
137 28915 : end subroutine abi_dhegv
138 : !!***
139 :
140 : !----------------------------------------------------------------------
141 :
142 : !!****f* m_abi_linalg/abi_chegv
143 : !! NAME
144 : !! abi_chegv
145 : !!
146 : !! FUNCTION
147 : !!
148 : !! INPUTS
149 : !!
150 : !! SOURCE
151 : !!
152 0 : subroutine abi_chegv(itype,jobz,uplo,n,a,lda,b,ldb,w)
153 :
154 : implicit none
155 :
156 : !Arguments ------------------------------------
157 : integer,intent(in) :: itype
158 : character(len=1), intent(in) :: jobz
159 : character(len=1), intent(in) :: uplo
160 : integer, intent(in) :: n,lda,ldb
161 : complex(sp), intent(inout) :: a(lda,*)
162 : complex(sp), intent(inout) :: b(ldb,*)
163 : real(sp), intent(out) :: w(n)
164 :
165 : !Local variables-------------------------------
166 : integer :: info,lwork
167 : real(sp),pointer :: rwork(:)
168 : complex(sp),pointer :: work(:)
169 : ! *********************************************************************
170 :
171 0 : ABI_CHECK(lapack_full_storage,"BUG(1) in abi_chegv (storage)!")
172 0 : ABI_CHECK(lapack_single_precision,"BUG(2) in abi_chegv (precision)!")
173 0 : ABI_CHECK(n<=eigen_c_maxsize,"BUG(3) in abi_chegv (maxsize)!")
174 :
175 0 : work => eigen_c_work ; rwork => eigen_c_rwork
176 0 : lwork=eigen_c_lwork
177 :
178 : !===== PLASMA
179 : !FDahm & LNGuyen (November 2012) :
180 : ! In Plasma v 2.4.6, eigen routines support only
181 : ! the eigenvalues computation (jobz=N) and not the
182 : ! full eigenvectors bases determination (jobz=V)
183 0 : if (ABI_LINALG_PLASMA_ISON.and.LSAME(jobz,'N')) then
184 : #if defined HAVE_LINALG_PLASMA
185 : if (eigen_c_lwork==0) then
186 : ABI_MALLOC(work,(n**2))
187 : end if
188 : call PLASMA_Alloc_Workspace_chegv(n,n,plasma_work,info)
189 : info = call PLASMA_chegv(itype,jobz_plasma(jobz),uplo_plasma(uplo),n,a,lda,b,ldb,w,&
190 : & plasma_work,c_loc(work),n)
191 : call PLASMA_Dealloc_handle(plasma_work,info)
192 : if (eigen_c_lwork==0) then
193 : ABI_FREE(work)
194 : end if
195 : #endif
196 :
197 : !===== LAPACK
198 : else
199 0 : if (eigen_c_lwork==0) then
200 0 : lwork=2*n-1
201 0 : ABI_MALLOC(work,(lwork))
202 : end if
203 0 : if (eigen_c_lrwork==0) then
204 0 : ABI_MALLOC(rwork,(3*n-2))
205 : end if
206 0 : call chegv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,info)
207 0 : if (eigen_c_lwork==0) then
208 0 : ABI_FREE(work)
209 : end if
210 0 : if (eigen_c_lrwork==0) then
211 0 : ABI_FREE(rwork)
212 : end if
213 : end if
214 :
215 0 : ABI_CHECK(info==0,"abi_chegv returned info!=0!")
216 :
217 0 : end subroutine abi_chegv
218 : !!***
219 :
220 : !----------------------------------------------------------------------
221 :
222 : !!****f* m_abi_linalg/abi_zhegv
223 : !! NAME
224 : !! abi_zhegv
225 : !!
226 : !! FUNCTION
227 : !!
228 : !! INPUTS
229 : !!
230 : !!
231 : !! SOURCE
232 :
233 0 : subroutine abi_zhegv(itype,jobz,uplo,n,a,lda,b,ldb,w)
234 :
235 : implicit none
236 :
237 : !Arguments ------------------------------------
238 : integer,intent(in) :: itype
239 : character(len=1), intent(in) :: jobz
240 : character(len=1), intent(in) :: uplo
241 : integer, intent(in) :: n,lda,ldb
242 : complex(dp), intent(inout) :: a(lda,*)
243 : complex(dp), intent(inout) :: b(ldb,*)
244 : real(dp), intent(out) :: w(n)
245 :
246 : !Local variables-------------------------------
247 : integer :: info,lwork
248 : real(dp),pointer :: rwork(:)
249 : complex(dp),pointer :: work(:)
250 : ! *********************************************************************
251 :
252 0 : ABI_CHECK(lapack_full_storage,"BUG(1) in abi_zhegv (storage)!")
253 0 : ABI_CHECK(lapack_double_precision,"BUG(2) in abi_zhegv (precision)!")
254 0 : ABI_CHECK(n<=eigen_z_maxsize,"BUG(3) in abi_zhegv (maxsize)!")
255 :
256 0 : work => eigen_z_work ; rwork => eigen_z_rwork
257 0 : lwork=eigen_z_lwork
258 :
259 : !===== PLASMA
260 : !FDahm & LNGuyen (November 2012) :
261 : ! In Plasma v 2.4.6, eigen routines support only
262 : ! the eigenvalues computation (jobz=N) and not the
263 : ! full eigenvectors bases determination (jobz=V)
264 0 : if (ABI_LINALG_PLASMA_ISON.and.LSAME(jobz,'N')) then
265 : #if defined HAVE_LINALG_PLASMA
266 : if (eigen_z_lwork==0) then
267 : ABI_MALLOC(work,(n**2))
268 : end if
269 : call PLASMA_Alloc_Workspace_zhegv(n,n,plasma_work,info)
270 : call PLASMA_zhegv(itype,jobz_plasma(jobz),uplo_plasma(uplo),n,a,lda,b,ldb,w,&
271 : & plasma_work,c_loc(work),n)
272 : call PLASMA_Dealloc_handle(plasma_work,info)
273 : if (eigen_z_lwork==0) then
274 : ABI_FREE(work)
275 : end if
276 : #endif
277 :
278 : !===== LAPACK
279 : else
280 0 : if (eigen_z_lwork==0) then
281 0 : lwork=2*n-1
282 0 : ABI_MALLOC(work,(lwork))
283 : end if
284 0 : if (eigen_z_lrwork==0) then
285 0 : ABI_MALLOC(rwork,(3*n-2))
286 : end if
287 0 : call zhegv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,info)
288 0 : if (eigen_z_lwork==0) then
289 0 : ABI_FREE(work)
290 : end if
291 0 : if (eigen_z_lrwork==0) then
292 0 : ABI_FREE(rwork)
293 : end if
294 : end if
295 :
296 0 : ABI_CHECK(info==0,"abi_zhegv returned info!=0!")
297 :
298 0 : end subroutine abi_zhegv
299 : !!***
300 :
301 : !----------------------------------------------------------------------
302 :
303 : !!****f* m_abi_linalg/abi_d2zhegvd
304 : !! NAME
305 : !! abi_d2zhegvd
306 : !!
307 : !! FUNCTION
308 : !! Generic divide-and-conquer generalized Hermitian eigensolver (HEGVD/SYGVD)
309 : !! with GPU support. Accepts real storage (real or complex-as-real via x_cplx).
310 : !! On GPU, calls abi_gpu_xhegvd_cptr. On CPU, self-manages work arrays.
311 : !!
312 : !! INPUTS
313 : !!
314 : !! SOURCE
315 :
316 0 : subroutine abi_d2zhegvd(itype, jobz, uplo, n, a, lda, b, ldb, w, info, x_cplx, gpu_option)
317 :
318 : !Arguments ------------------------------------
319 : integer, intent(in) :: itype
320 : character(len=1), intent(in) :: jobz
321 : character(len=1), intent(in) :: uplo
322 : integer, intent(in) :: n, lda, ldb
323 : real(dp), target, intent(inout) :: a(:,:)
324 : real(dp), target, intent(inout) :: b(:,:)
325 : real(dp), target, intent(out) :: w(:,:)
326 : integer, intent(out) :: info
327 : !Optionals -----------------------------------
328 : integer, intent(in), optional :: x_cplx
329 : integer, intent(in), optional :: gpu_option
330 :
331 : !Local variables-------------------------------
332 : integer :: cplx_, gpu_option_
333 : integer :: lwork, lrwork, liwork
334 0 : real(dp), pointer :: rwork(:)
335 0 : complex(dp), pointer :: cwork(:)
336 0 : integer, pointer :: iwork(:)
337 : real(dp) :: rwork_query(1)
338 : complex(dp) :: cwork_query(1)
339 : integer :: iwork_query(1)
340 :
341 : ! *********************************************************************
342 :
343 0 : cplx_=1 ; if(PRESENT(x_cplx)) cplx_ = x_cplx
344 0 : gpu_option_=ABI_GPU_DISABLED ; if(PRESENT(gpu_option)) gpu_option_ = gpu_option
345 :
346 : #if defined(DEBUG_VERBOSE) && defined(HAVE_OPENMP_OFFLOAD)
347 : if ( gpu_option_ == ABI_GPU_OPENMP ) then
348 : ABI_CHECK(xomp_target_is_present(c_loc(a)), "Array isn't mapped on GPU")
349 : ABI_CHECK(xomp_target_is_present(c_loc(b)), "Array isn't mapped on GPU")
350 : ABI_CHECK(xomp_target_is_present(c_loc(w)), "Array isn't mapped on GPU")
351 : end if
352 : #endif
353 :
354 0 : if(gpu_option_/=ABI_GPU_DISABLED) then
355 0 : if(gpu_option_==ABI_GPU_OPENMP) then
356 : #ifdef HAVE_OPENMP_OFFLOAD
357 : !$OMP TARGET DATA USE_DEVICE_ADDR(a,b,w)
358 : call abi_gpu_xhegvd_cptr(cplx_, itype, jobz, uplo, n, c_loc(a), lda, c_loc(b), ldb, c_loc(w), info)
359 : !$OMP END TARGET DATA
360 : #endif
361 : else
362 0 : call abi_gpu_xhegvd_cptr(cplx_, itype, jobz, uplo, n, c_loc(a), lda, c_loc(b), ldb, c_loc(w), info)
363 : end if
364 : else
365 0 : if(cplx_ == 2) then
366 0 : lwork=-1 ; lrwork=-1 ; liwork=-1
367 0 : call zhegvd(itype, jobz, uplo, n, a, lda, b, ldb, w, cwork_query, lwork, rwork_query, lrwork, iwork_query, liwork, info)
368 0 : lwork=int(cwork_query(1)) ; lrwork=int(rwork_query(1)) ; liwork=iwork_query(1)
369 0 : ABI_MALLOC(cwork, (lwork))
370 0 : ABI_MALLOC(rwork, (lrwork))
371 0 : ABI_MALLOC(iwork, (liwork))
372 0 : call zhegvd(itype, jobz, uplo, n, a, lda, b, ldb, w, cwork, lwork, rwork, lrwork, iwork, liwork, info)
373 0 : ABI_FREE(cwork) ; ABI_FREE(rwork) ; ABI_FREE(iwork)
374 : else
375 0 : lwork=-1 ; liwork=-1
376 0 : call dsygvd(itype, jobz, uplo, n, a, lda, b, ldb, w, rwork_query, lwork, iwork_query, liwork, info)
377 0 : lwork=int(rwork_query(1)) ; liwork=iwork_query(1)
378 0 : ABI_MALLOC(rwork, (lwork))
379 0 : ABI_MALLOC(iwork, (liwork))
380 0 : call dsygvd(itype, jobz, uplo, n, a, lda, b, ldb, w, rwork, lwork, iwork, liwork, info)
381 0 : ABI_FREE(rwork) ; ABI_FREE(iwork)
382 : end if
383 : end if
384 :
385 0 : ABI_CHECK(info==0,"abi_d2zhegvd returned info!=0!")
386 :
387 0 : end subroutine abi_d2zhegvd
388 : !!***
389 :
390 : !----------------------------------------------------------------------
391 :
392 : !!****f* m_abi_linalg/abi_zhegvd_2d
393 : !! NAME
394 : !! abi_zhegvd_2d
395 : !!
396 : !! FUNCTION
397 : !! Divide-and-conquer generalized complex Hermitian eigensolver (ZHEGVD) with
398 : !! GPU support. Accepts complex(dp) 2D matrices. On GPU, calls
399 : !! abi_gpu_xhegvd_cptr. On CPU, self-manages work arrays.
400 : !!
401 : !! INPUTS
402 : !!
403 : !! SOURCE
404 :
405 0 : subroutine abi_zhegvd_2d(itype, jobz, uplo, n, a, lda, b, ldb, w, info, gpu_option)
406 :
407 : !Arguments ------------------------------------
408 : integer, intent(in) :: itype
409 : character(len=1), intent(in) :: jobz
410 : character(len=1), intent(in) :: uplo
411 : integer, intent(in) :: n, lda, ldb
412 : complex(dp), target, intent(inout) :: a(:,:)
413 : complex(dp), target, intent(inout) :: b(:,:)
414 : real(dp), target, intent(out) :: w(:,:)
415 : integer, intent(out) :: info
416 : !Optionals -----------------------------------
417 : integer, intent(in), optional :: gpu_option
418 :
419 : !Local variables-------------------------------
420 : integer :: gpu_option_
421 : integer :: lwork, lrwork, liwork
422 0 : complex(dp), pointer :: cwork(:)
423 0 : real(dp), pointer :: rwork(:)
424 0 : integer, pointer :: iwork(:)
425 : complex(dp) :: cwork_query(1)
426 : real(dp) :: rwork_query(1)
427 : integer :: iwork_query(1)
428 :
429 : ! *********************************************************************
430 :
431 0 : gpu_option_=ABI_GPU_DISABLED ; if(PRESENT(gpu_option)) gpu_option_ = gpu_option
432 :
433 : #if defined(DEBUG_VERBOSE) && defined(HAVE_OPENMP_OFFLOAD)
434 : if ( gpu_option_ == ABI_GPU_OPENMP ) then
435 : ABI_CHECK(xomp_target_is_present(c_loc(a)), "Array isn't mapped on GPU")
436 : ABI_CHECK(xomp_target_is_present(c_loc(b)), "Array isn't mapped on GPU")
437 : ABI_CHECK(xomp_target_is_present(c_loc(w)), "Array isn't mapped on GPU")
438 : end if
439 : #endif
440 :
441 0 : if(gpu_option_/=ABI_GPU_DISABLED) then
442 0 : if(gpu_option_==ABI_GPU_OPENMP) then
443 : #ifdef HAVE_OPENMP_OFFLOAD
444 : !$OMP TARGET DATA USE_DEVICE_ADDR(a,b,w)
445 : call abi_gpu_xhegvd_cptr(2, itype, jobz, uplo, n, c_loc(a), lda, c_loc(b), ldb, c_loc(w), info)
446 : !$OMP END TARGET DATA
447 : #endif
448 : else
449 0 : call abi_gpu_xhegvd_cptr(2, itype, jobz, uplo, n, c_loc(a), lda, c_loc(b), ldb, c_loc(w), info)
450 : end if
451 : else
452 0 : lwork=-1 ; lrwork=-1 ; liwork=-1
453 0 : call zhegvd(itype, jobz, uplo, n, a, lda, b, ldb, w, cwork_query, lwork, rwork_query, lrwork, iwork_query, liwork, info)
454 0 : lwork=int(cwork_query(1)) ; lrwork=int(rwork_query(1)) ; liwork=iwork_query(1)
455 0 : ABI_MALLOC(cwork, (lwork))
456 0 : ABI_MALLOC(rwork, (lrwork))
457 0 : ABI_MALLOC(iwork, (liwork))
458 0 : call zhegvd(itype, jobz, uplo, n, a, lda, b, ldb, w, cwork, lwork, rwork, lrwork, iwork, liwork, info)
459 0 : ABI_FREE(cwork) ; ABI_FREE(rwork) ; ABI_FREE(iwork)
460 : end if
461 :
462 0 : ABI_CHECK(info==0,"abi_zhegvd_2d returned info!=0!")
463 :
464 0 : end subroutine abi_zhegvd_2d
465 : !!***
|
