Line data Source code
1 : !!****f* ABINIT/m_trace_estimation
2 : !! NAME
3 : !! m_trace_estimation
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines used to compute the stochastic Lanczos trace estimation.
7 : !! It also computes cummulative eigenvalue counts using differences of trace estimation
8 : !! on consecutive intervals.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2018-2026 ABINIT group (IML)
12 : !! This file is distributed under the terms of the
13 : !! gnu general public license, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : ! nvtx related macro definition
26 : #include "nvtx_macros.h"
27 :
28 : module m_trace_estimation
29 :
30 : use defs_basis
31 : use defs_abitypes
32 : use m_abicore
33 : use m_errors
34 : use m_time, only : timab
35 : use m_sort, only: sort_dp
36 :
37 : use m_cgtools
38 : use m_xg
39 : use m_xgTransposer
40 :
41 : use m_chebfi2
42 : use m_polynomial_filter
43 : use m_slice_task, only: matrixInfo_t
44 :
45 : use m_xmpi
46 : use m_xomp
47 : #ifdef HAVE_OPENMP
48 : use omp_lib
49 : #endif
50 :
51 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_YAKL)
52 : use m_gpu_toolbox, only : CPU_DEVICE_ID, gpu_device_synchronize
53 : #endif
54 :
55 : #if defined(HAVE_GPU_MARKERS)
56 : use m_nvtx_data
57 : #endif
58 :
59 : implicit none
60 :
61 : private
62 :
63 : ! Timers
64 : !---------------------------------------------------
65 : integer, parameter :: tim_swap = 1761
66 : integer, parameter :: tim_RR_q = 1759
67 : integer, parameter :: tim_barrier = 1764
68 : integer, parameter :: tim_copy = 1765
69 : integer, parameter :: tim_getAX_BX = 1754
70 : integer, parameter :: tim_invovl = 1755
71 : integer, parameter :: tim_lanczos = 2167
72 : integer, parameter :: tim_trace = 2168
73 :
74 : ! Public methods
75 : !-------------------------------------------------
76 : public :: computeBLanczos ! for lower bound estimation
77 : public :: get_eigenvalue_count ! count eigenvalues for lowpass intervals
78 : public :: computeTraceEstimation ! compute moments etc
79 : public :: smallestTridiagEigenpair ! used in slice_cprj (experimental)
80 :
81 : CONTAINS
82 : !=====================================================================
83 : !!***
84 :
85 : !!****f* m_trace_estimation/computeBLanczos
86 : !! NAME
87 : !! computeBLanczos
88 : !!
89 : !! FUNCTION
90 : !! B-Lanczos three-term recurrence (using B-inner product).
91 : !! Performs k Lanczos iterations on a column vector.
92 : !!
93 : !! SOURCE
94 :
95 18 : subroutine computeBLanczos(minfo, paw, getAX_BX, getBm1X, k, lambda_min, res_norm)
96 :
97 : implicit none
98 :
99 : type(matrixInfo_t), intent(in) :: minfo
100 : logical, intent(in) :: paw
101 : integer, intent(in) :: k
102 : real(dp), intent(out) :: lambda_min, res_norm
103 : interface
104 : subroutine getAX_BX(X,AX,BX)
105 : use m_xg, only : xgBlock_t
106 : type(xgBlock_t), intent(inout) :: X
107 : type(xgBlock_t), intent(inout) :: AX
108 : type(xgBlock_t), intent(inout) :: BX
109 : end subroutine getAX_BX
110 : end interface
111 : interface
112 : subroutine getBm1X(X,Bm1X)
113 : use m_xg, only : xgBlock_t
114 : type(xgBlock_t), intent(inout) :: X
115 : type(xgBlock_t), intent(inout) :: Bm1X
116 : end subroutine getBm1X
117 : end interface
118 :
119 : type(xg_t) :: W_vcol
120 : type(xg_t) :: W_dot
121 : type(xgBlock_t) :: q, v, Bv, Bm1v, qprev
122 : type(xgBlock_t) :: dot_qTBv, dot_qTv, dot_vTBv
123 36 : real(dp) :: alpha(k), beta(k-1)
124 : real(dp) :: beta_prev
125 : real(dp) :: norml_q
126 : real(dp), pointer :: dot_qTBv_layout(:,:) => null()
127 : real(dp), pointer :: dot_qTv_layout(:,:) => null()
128 : real(dp), pointer :: dot_vTBv_layout(:,:) => null()
129 18 : real(dp), allocatable :: v_min(:)
130 :
131 : integer :: j
132 : integer :: rank
133 : integer :: me_g0, space, tot_spacedim, gpu_option
134 : real(dp) :: tsec(2)
135 :
136 : ! *********************************************************************
137 :
138 18 : call timab(tim_lanczos,1,tsec)
139 :
140 18 : tot_spacedim = minfo%total_spacedim
141 18 : gpu_option = minfo%gpu_option
142 18 : space = minfo%space
143 18 : me_g0 = minfo%me_g0
144 18 : if (minfo%paral_kgb==1) then
145 16 : me_g0 = minfo%me_g0_fft
146 : end if
147 18 : rank = xmpi_comm_rank(minfo%spacecom)
148 18 : beta_prev = 0.0_dp
149 :
150 54 : ABI_MALLOC(v_min, (tot_spacedim))
151 :
152 18 : write(std_out,*) 'Lanczos in rank=', rank; flush(std_out)
153 :
154 : ! workspace size (npw,5)
155 18 : call xg_init(W_vcol, space, tot_spacedim, 5, xmpi_comm_null, me_g0=me_g0, gpu_option=gpu_option)
156 18 : call xgBlock_setBlock(W_vcol%self, q, tot_spacedim, 1) ! q
157 18 : call xgBlock_setBlock(W_vcol%self, v, tot_spacedim, 1, fcol=2) ! Aq
158 18 : call xgBlock_setBlock(W_vcol%self, Bv, tot_spacedim, 1, fcol=3) ! Bq
159 18 : call xgBlock_setBlock(W_vcol%self, Bm1v, tot_spacedim, 1, fcol=4) ! Bm1 v
160 18 : call xgBlock_setBlock(W_vcol%self, qprev, tot_spacedim, 1, fcol=5) ! q_prev
161 :
162 18 : call xg_init(W_dot, space, 1, 3, xmpi_comm_null, me_g0=me_g0, gpu_option=gpu_option)
163 18 : call xgBlock_setBlock(W_dot%self, dot_qTBv, 1, 1)
164 18 : call xgBlock_setBlock(W_dot%self, dot_qTv, 1, 1, fcol=2)
165 18 : call xgBlock_setBlock(W_dot%self, dot_vTBv, 1, 1, fcol=3)
166 :
167 : ! q = random column vector
168 18 : call xgBlock_colwiseRandom(q, rank, 1)
169 : !write(std_out,*) 'Random id=', xgBlock_getid(q)
170 : !flush(std_out)
171 :
172 : ! Bv = B * q / norml_q
173 18 : call timab(tim_getAX_BX,1,tsec)
174 : ABI_NVTX_START_RANGE(NVTX_SLICE_GET_AX_BX)
175 18 : call getAX_BX(q, v, Bv)
176 18 : call xgBlock_zero_im_g0(v) ! v stores Aq
177 18 : call xgBlock_zero_im_g0(Bv) ! Bv stores Bq
178 : ABI_NVTX_END_RANGE()
179 18 : call timab(tim_getAX_BX,2,tsec)
180 :
181 18 : call xgBlock_colwiseDotProduct(q, Bv, dot_qTBv)
182 18 : call xgBlock_reverseMap(dot_qTBv,dot_qTBv_layout,rows=1,cols=1)
183 :
184 18 : norml_q = 1.d0 / sqrt(dot_qTBv_layout(1,1))
185 18 : call xgBlock_scale(q, norml_q, 1)
186 18 : call xgBlock_scale(v, norml_q, 1)
187 18 : call xgBlock_scale(Bv, norml_q, 1)
188 :
189 558 : alpha = 0.d0
190 540 : beta = 0.d0
191 :
192 558 : do j = 1, k
193 :
194 540 : if (j>1) then ! rewrite v
195 : ! v = A * q
196 522 : call timab(tim_getAX_BX,1,tsec)
197 : ABI_NVTX_START_RANGE(NVTX_SLICE_GET_AX_BX)
198 522 : call getAX_BX(q, v, Bv)
199 522 : call xgBlock_zero_im_g0(v) ! v stores Aq
200 522 : call xgBlock_zero_im_g0(Bv) ! Bv stores Bq
201 : ABI_NVTX_END_RANGE()
202 522 : call timab(tim_getAX_BX,2,tsec)
203 : end if
204 :
205 : ! alpha_j = q^T * Aq
206 540 : call xgBlock_colwiseDotProduct(q, v, dot_qTv)
207 540 : call xgBlock_reverseMap(dot_qTv,dot_qTv_layout,rows=1,cols=1)
208 540 : alpha(j) = dot_qTv_layout(1,1)
209 :
210 : ! v = B^{-1} * A * q
211 540 : if (paw) then
212 540 : call timab(tim_invovl, 1, tsec)
213 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_BM1X)
214 540 : call getBm1X(v, Bm1v)
215 : ABI_NVTX_END_RANGE()
216 540 : call timab(tim_invovl, 2, tsec)
217 540 : call timab(tim_copy, 1, tsec)
218 540 : call xgBlock_copy(Bm1v, v)
219 540 : call timab(tim_copy, 2, tsec)
220 : end if
221 :
222 : ! v = B^{-1} A q - alpha q - beta_prev q_prev
223 540 : call xgBlock_saxpy(v, -1.d0 * alpha(j), q)
224 540 : if (j > 1) then
225 522 : call xgBlock_saxpy(v, -beta_prev, qprev)
226 : end if
227 :
228 : ! Compute beta_j if j<k
229 558 : if (j < k) then
230 :
231 : ! Bv = B * v
232 522 : call timab(tim_getAX_BX,1,tsec)
233 : ABI_NVTX_START_RANGE(NVTX_SLICE_GET_AX_BX)
234 522 : call getAX_BX(v, Bm1v, Bv)
235 522 : call xgBlock_zero_im_g0(Bm1v) ! Bm1v dummy workspace
236 522 : call xgBlock_zero_im_g0(Bv) ! Bv
237 : ABI_NVTX_END_RANGE()
238 522 : call timab(tim_getAX_BX,2,tsec)
239 :
240 522 : call xgBlock_colwiseDotProduct(v, Bv, dot_vTBv)
241 522 : call xgBlock_reverseMap(dot_vTBv,dot_vTBv_layout,rows=1,cols=1)
242 522 : beta(j) = sqrt(dot_vTBv_layout(1,1))
243 :
244 : ! Update q_prev, q=v/beta_j, beta_prev
245 522 : call xgBlock_copy(q, qprev)
246 522 : call xgBlock_scale(v, 1.d0/beta(j), 1)
247 522 : call xgBlock_copy(v, q)
248 522 : beta_prev = beta(j)
249 : end if
250 : end do
251 :
252 : ! Diagonalize T (always on CPU)
253 18 : call smallestTridiagEigenpair(k, alpha, beta, lambda_min, v_min)
254 :
255 : ! residual norm using Lanczos shortcut
256 18 : res_norm = abs(beta(k-1)*v_min(k))
257 :
258 18 : call xg_free(W_vcol)
259 18 : call xg_free(W_dot)
260 18 : ABI_FREE(v_min)
261 :
262 18 : call timab(tim_lanczos,2,tsec)
263 :
264 36 : end subroutine computeBLanczos
265 : !!***
266 :
267 : !----------------------------------------------------------------------
268 :
269 : !!****f* m_trace_estimation/computeChebyshevMoments
270 : !! NAME
271 : !! computeChebyshevMoments
272 : !!
273 : !! FUNCTION
274 : !! Compute Chebyshev moments up to maximal degree all centered in [A,B)
275 : !! Upper bound is computed as Rayleigh quotient
276 : !!
277 : !! OUTPUT
278 : !! M_n = <X, f_n(B^{-1}AX) X> for n=1,..,ndeg_filter_max
279 : !!
280 : !! SOURCE
281 :
282 18 : subroutine computeChebyshevMoments(minfo, tolerance, ecut, paw, X0, getAX_BX, getBm1X, &
283 18 : lambda_minus, lambda_plus, ndeg_filter, cheby_moments)
284 :
285 : implicit none
286 :
287 : type(matrixInfo_t), intent(in) :: minfo
288 : real(dp), intent(in) :: ecut, tolerance
289 : logical, intent(in) :: paw
290 : type(xgBlock_t), intent(inout) :: X0
291 : integer, intent(in) :: ndeg_filter
292 : real(dp), intent(in) :: lambda_minus, lambda_plus
293 : complex(dp), intent(out) :: cheby_moments(:,:)
294 : interface
295 : subroutine getAX_BX(X,AX,BX)
296 : use m_xg, only : xgBlock_t
297 : type(xgBlock_t), intent(inout) :: X
298 : type(xgBlock_t), intent(inout) :: AX
299 : type(xgBlock_t), intent(inout) :: BX
300 : end subroutine getAX_BX
301 : end interface
302 : interface
303 : subroutine getBm1X(X,Bm1X)
304 : use m_xg, only : xgBlock_t
305 : type(xgBlock_t), intent(inout) :: X
306 : type(xgBlock_t), intent(inout) :: Bm1X
307 : end subroutine getBm1X
308 : end interface
309 :
310 : integer :: ideg
311 : integer :: nband, tot_spacedim, space, me_g0, gpu_option
312 : real(dp) :: center, radius
313 : real(dp) :: one_over_r
314 : real(dp) :: two_over_r
315 : real(dp) :: tsec(2)
316 : type(xg_t) :: Moments
317 : type(xg_t) :: X0_backup
318 : type(xgBlock_t) :: Moment_ideg
319 18 : type(chebfi_t) :: chebfi
320 : complex(dp), pointer :: momvals(:,:) => null()
321 :
322 : ! *********************************************************************
323 :
324 : ! todo add timer
325 : ! tim_cheby_moments
326 :
327 18 : gpu_option = minfo%gpu_option
328 18 : tot_spacedim = minfo%total_spacedim
329 18 : space = minfo%space
330 18 : me_g0 = minfo%me_g0
331 18 : if (minfo%paral_kgb==1) then
332 16 : me_g0 = minfo%me_g0_fft
333 : end if
334 18 : nband = cols(X0)
335 :
336 : ! Moment workspace size (1, ndeg+1)
337 18 : call xg_init(Moments, space, 1, ndeg_filter+1, gpu_option=gpu_option) ! M_n=<X0,f_n(A)X0>
338 18 : call xg_init(X0_backup, space, tot_spacedim, nband, minfo%spacecom, me_g0=me_g0, gpu_option=gpu_option) ! X0
339 :
340 : ! Initialize chebfi object in MPI Colsrows distribution
341 : call chebfi_init(chebfi,nband,tot_spacedim,tolerance,ecut,minfo%paral_kgb,&
342 : nband,ndeg_filter,0,minfo%space,1,xmpi_comm_null,minfo%me_g0,minfo%me_g0_fft,&
343 : paw,minfo%comm_rows,minfo%comm_cols,0,1.d0,0.d0,minfo%gpu_option,&
344 : gpu_kokkos_nthrd=minfo%gpu_kokkos_nthrd,gpu_thread_limit=minfo%gpu_thread_limit,&
345 18 : from_linalg=.false.)
346 :
347 : ! Initialize Chebyshev recursion
348 18 : call xgBlock_copy(X0, X0_backup%self)
349 18 : chebfi%xXColsRows = X0
350 :
351 : ! Compute A*Psi
352 18 : call timab(tim_getAX_BX,1,tsec)
353 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_GET_AX_BX)
354 18 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
355 18 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
356 18 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
357 : ABI_NVTX_END_RANGE()
358 18 : call timab(tim_getAX_BX,2,tsec)
359 :
360 : ! Initialize Chebyshev moment at k=1
361 18 : call xgBlock_setBlock(Moments%self, Moment_ideg, 1, 1, fcol=1)
362 18 : call xgBlock_dot(X0_backup%self, chebfi%xXColsRows, Moment_ideg)
363 :
364 : ! Spectral interval to be amplified scaled to [-1,1)
365 18 : center = (lambda_plus + lambda_minus)/2.d0
366 18 : radius = (lambda_plus - lambda_minus)/2.d0
367 18 : one_over_r = 1.d0/radius
368 18 : two_over_r = 2.d0/radius
369 :
370 198 : do ideg = 0, ndeg_filter - 1
371 :
372 : !chebfi%paw = .false.
373 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_NEXT_ORDER)
374 180 : call chebfi_computeNextOrderChebfiPolynom(chebfi, ideg, center, one_over_r, two_over_r, getBm1X)
375 : ABI_NVTX_END_RANGE()
376 : !chebfi%paw = .true.
377 :
378 : ! chebfi%xXColsRows = f_ideg X0
379 : ABI_NVTX_START_RANGE(NVTX_CHEBFI2_SWAP_BUF)
380 180 : call timab(tim_swap,1,tsec)
381 180 : call chebfi_swapInnerBuffers(chebfi, tot_spacedim, nband)
382 180 : call timab(tim_swap,2,tsec)
383 : ABI_NVTX_END_RANGE()
384 :
385 : !A * Psi
386 180 : call timab(tim_getAX_BX,1,tsec)
387 : ABI_NVTX_START_RANGE(NVTX_SLICE_GET_AX_BX)
388 180 : call getAX_BX(chebfi%xXColsRows, chebfi%xAXColsRows, chebfi%xBXColsRows)
389 180 : call xgBlock_zero_im_g0(chebfi%xAXColsRows)
390 180 : call xgBlock_zero_im_g0(chebfi%xBXColsRows)
391 : ABI_NVTX_END_RANGE()
392 180 : call timab(tim_getAX_BX,2,tsec)
393 :
394 : ! M_ideg = < X0, f_ideg X0 >_B
395 180 : call xgBlock_setBlock(Moments%self, Moment_ideg, 1, 1, fcol=ideg+2)
396 198 : call xgBlock_dot(X0_backup%self, chebfi%xXColsRows, Moment_ideg)
397 :
398 : end do
399 :
400 18 : call xgBlock_reverseMap(Moments%self, momvals, 1, ndeg_filter+1)
401 414 : cheby_moments(:,:) = momvals(:,:)
402 :
403 : ! Free memory
404 18 : call xg_free(Moments)
405 18 : call chebfi_free(chebfi)
406 18 : call xg_free(X0_backup)
407 :
408 36 : end subroutine computeChebyshevMoments
409 : !!***
410 :
411 : !----------------------------------------------------------------------
412 :
413 : !!****f* m_trace_estimation/computeTraceEstimation
414 : !! NAME
415 : !! computeTraceEstimation
416 : !!
417 : !! FUNCTION
418 : !! Compute Girard-Hutchinson trace estimator
419 : !! ndeg_filter -> number of Chebyshev moments
420 : !! m_probe -> number of stochastic probes
421 : !! moments(1:ndeg_filter) -> Chebyshev moments computed with probes
422 : !! and B-inner product
423 : !!
424 : !! SOURCE
425 :
426 18 : subroutine computeTraceEstimation(minfo, ecut, paw, tolerance, getAX_BX, getBm1X, ndeg_filter, &
427 18 : m_probe, min_low_bound, moments)
428 :
429 : implicit none
430 :
431 : type(matrixInfo_t), intent(in) :: minfo
432 : integer, intent(in) :: ndeg_filter, m_probe
433 : logical, intent(in) :: paw
434 : real(dp), intent(in) :: ecut, min_low_bound, tolerance
435 : real(dp), intent(inout) :: moments(:)
436 : interface
437 : subroutine getAX_BX(X,AX,BX)
438 : use m_xg, only : xgBlock_t
439 : type(xgBlock_t), intent(inout) :: X
440 : type(xgBlock_t), intent(inout) :: AX
441 : type(xgBlock_t), intent(inout) :: BX
442 : end subroutine getAX_BX
443 : end interface
444 : interface
445 : subroutine getBm1X(X,Bm1X)
446 : use m_xg, only : xgBlock_t
447 : type(xgBlock_t), intent(inout) :: X
448 : type(xgBlock_t), intent(inout) :: Bm1X
449 : end subroutine getBm1X
450 : end interface
451 :
452 : !integer :: nband
453 : integer :: tot_spacedim, gpu_option
454 : integer :: ierr
455 : integer :: me_g0
456 : integer :: jcol
457 : integer :: my_rank
458 : integer :: seed
459 : integer :: m_probe_tot
460 : integer :: num_moments
461 : type(xg_t) :: X_probe
462 : complex(dp), allocatable :: cheby_moments(:,:)
463 : real(dp) :: tsec(2)
464 :
465 : ! *********************************************************************
466 :
467 18 : call timab(tim_trace,1,tsec)
468 :
469 18 : tot_spacedim = minfo%total_spacedim
470 18 : gpu_option = minfo%gpu_option
471 18 : my_rank = xmpi_comm_rank(minfo%spacecom)
472 18 : me_g0 = minfo%me_g0
473 18 : if (minfo%paral_kgb==1) then
474 16 : me_g0 = minfo%me_g0_fft
475 : end if
476 :
477 18 : num_moments = ndeg_filter + 1
478 54 : ABI_MALLOC(cheby_moments, (1, num_moments) )
479 :
480 : ! total number of probes is m_probes * number of MPI processes
481 : call xg_init(X_probe, minfo%space, tot_spacedim, m_probe, minfo%spacecom, &
482 18 : me_g0=me_g0, gpu_option=gpu_option)
483 :
484 : ! Define random isotropic probes (unit variance NOT unit norm!)
485 108 : do jcol = 1, m_probe
486 90 : seed = my_rank*(m_probe+10)+jcol ! seed depends on column index
487 108 : call xgBlock_colwiseRandomRademacher(X_probe%self, seed, jcol)
488 : end do
489 :
490 18 : write(std_out,*) 'moments are scaled in'
491 18 : write(std_out,*) min_low_bound, ecut
492 18 : flush(std_out)
493 :
494 : ! Lowpass scan - coarse resolution with low degree
495 : ! Chebyshev moments in maximal [a,ecut)
496 : ! todo <X_probe, f(A) X_probe> (trace) and <X0, f(A) X_probe> (principal angles)
497 : ! maybe <xX, f(A)X_probe> is useful for principal angles and column selection
498 : ! in that case incorporate it in the loop
499 : call computeChebyshevMoments(minfo, tolerance, ecut, paw, X_probe%self, &
500 18 : getAX_BX, getBm1X, min_low_bound, ecut, ndeg_filter, cheby_moments)
501 :
502 : ! Sum real part of moments and divide by number of probes
503 : ! moments(k) = 1/Nv * Sum_{i=1}^Nv v_i^T T_k(A)v_i
504 18 : m_probe_tot = m_probe
505 18 : call xmpi_sum(m_probe_tot, minfo%spacecom, ierr)
506 : !moments(1:num_moments) = (/ (real(cheby_moments(1,k)), k=1,num_moments) /)
507 216 : moments = real(cheby_moments(1,:))
508 18 : call xmpi_sum(moments, minfo%spacecom, ierr)
509 216 : moments(1:num_moments) = moments(1:num_moments)/m_probe_tot
510 :
511 : !write(std_out,*) 'moments k=0=', real(cheby_moments(1,1))
512 : !write(std_out,*) 'moments k=1=', real(cheby_moments(1,2))
513 : !write(std_out,*) 'moments k=3=', real(cheby_moments(1,3))
514 : !flush(std_out)
515 :
516 : ! Free memory
517 18 : ABI_FREE(cheby_moments)
518 18 : call xg_free(X_probe)
519 :
520 18 : call timab(tim_trace,2,tsec)
521 :
522 36 : end subroutine computeTraceEstimation
523 : !!***
524 :
525 : !----------------------------------------------------------------------
526 :
527 : !!****f* m_trace_estimation/get_eigenvalue_count
528 : !! NAME
529 : !! get_eigenvalue_count
530 : !!
531 : !! SOURCE
532 :
533 612 : function get_eigenvalue_count(b, moments, work) result(mass)
534 :
535 : implicit none
536 :
537 : real(dp), intent(in) :: b
538 : real(dp), intent(in) :: moments(:)
539 : real(dp), intent(inout) :: work(:)
540 : real(dp) :: mass
541 : integer :: ndeg_filter
542 : !real(dp) :: sigma, alpha
543 : !integer :: Ngrid
544 :
545 612 : ndeg_filter = size(moments)-1
546 :
547 : ! Erf damping coefficients
548 : !sigma = 4.d0 / ndeg_filter
549 : !Ngrid = 500
550 : !work = erf_step_coeffs(b, ndeg_filter, sigma, Ngrid)
551 :
552 : !alpha = 20
553 : !Ngrid = 500
554 : !work = smooth_step_coeffs(b, ndeg_filter, alpha, Ngrid)
555 :
556 : ! Steep Lanczos
557 7344 : work = lanczos_step_coeffs(b, ndeg_filter)
558 :
559 7344 : mass = dot_product(work, moments)
560 :
561 612 : end function get_eigenvalue_count
562 : !!***
563 :
564 : !----------------------------------------------------------------------
565 :
566 : !!****f* m_trace_estimation/smallestTridiagEigenpair
567 : !! NAME
568 : !! smallestTridiagEigenpair
569 : !!
570 : !! FUNCTION
571 : !! Smallest eigenvalue of symmetric tridiagonal
572 : !!
573 : !! SOURCE
574 :
575 18 : subroutine smallestTridiagEigenpair(n, d, e, lambda, v)
576 : implicit none
577 : integer, intent(in) :: n
578 : real(dp), intent(in) :: d(n), e(n-1)
579 : real(dp), intent(out) :: lambda
580 : real(dp), intent(out) :: v(n)
581 :
582 : ! Local copies (DSTEVX overwrites input)
583 36 : real(dp) :: dloc(n), eloc(n-1)
584 18 : real(dp), allocatable :: z(:,:), work(:)
585 18 : integer, allocatable :: iwork(:), ifail(:)
586 : integer :: info, m
587 :
588 : ! *********************************************************************
589 :
590 558 : dloc = d
591 540 : eloc = e
592 :
593 54 : ABI_MALLOC(z, (n,1))
594 54 : ABI_MALLOC(work, (5*n))
595 54 : ABI_MALLOC(iwork, (5*n))
596 54 : ABI_MALLOC(ifail, (n))
597 :
598 : ! DSTEVX computes selected eigenpairs (here smallest: index 1)
599 : call dstevx('V', 'I', n, dloc, eloc, 0.0d0, 0.0d0, 1, 1, 1.0d-12, m, dloc, z, &
600 18 : n, work, iwork, ifail, info)
601 :
602 18 : if (info /= 0) then
603 0 : ABI_ERROR('DSTEVX failed')
604 : end if
605 :
606 18 : lambda = dloc(1)
607 558 : v = z(:,1)
608 :
609 18 : ABI_FREE(z)
610 18 : ABI_FREE(work)
611 18 : ABI_FREE(iwork)
612 18 : ABI_FREE(ifail)
613 :
614 18 : end subroutine smallestTridiagEigenpair
615 : !!***
616 :
617 : end module m_trace_estimation
618 : !!***
|
