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