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