Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwls_utility
2 : : !! NAME
3 : : !! m_gwls_utility
4 : : !!
5 : : !! FUNCTION
6 : : !! .
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
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 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : :
23 : : !---------------------------------------------------------------------
24 : : ! Utility modules, which are not directly related to the
25 : : ! problem to be solved.
26 : : !---------------------------------------------------------------------
27 : :
28 : : module m_gwls_utility
29 : : !----------------------------------------------------------------------------------------------------
30 : : ! This module contains useful functions and subroutines, which are not necessarily numerical
31 : : ! in nature, for timing, opening files, etc...
32 : : !----------------------------------------------------------------------------------------------------
33 : :
34 : : ! abinit modules
35 : : use defs_basis
36 : : use m_abicore
37 : : use m_xmpi
38 : :
39 : : use m_abi_linalg
40 : : use m_io_tools, only : get_unit
41 : :
42 : : implicit none
43 : :
44 : : private
45 : :
46 : : complex(dp), public, parameter :: cmplx_i = (0.0_dp,1.0_dp)
47 : : complex(dp), public, parameter :: cmplx_1 = (1.0_dp,0.0_dp)
48 : : complex(dp), public, parameter :: cmplx_0 = (0.0_dp,0.0_dp)
49 : :
50 : : logical, public :: master_debug
51 : : character(len=100), public :: files_status_new='new'
52 : : character(len=100), public :: files_status_old='unknown'
53 : :
54 : : public :: driver_invert_positive_definite_hermitian_matrix
55 : : public :: ritz_analysis_general, orthogonalize
56 : : public :: complex_vector_product
57 : : !!***
58 : : contains
59 : :
60 : : !!****f* m_gwls_utility/complex_vector_product
61 : : !! NAME
62 : : !! complex_vector_product
63 : : !!
64 : : !! FUNCTION
65 : : !! .
66 : : !!
67 : : !! INPUTS
68 : : !!
69 : : !! OUTPUT
70 : : !!
71 : : !!
72 : : !!
73 : : !! SOURCE
74 : :
75 : 9848 : complex(dp) function complex_vector_product(v1,v2,l)
76 : : !--------------------------------------------------------------------------
77 : : ! This function computes the vector product of two complex vectors.
78 : : !--------------------------------------------------------------------------
79 : : integer, intent(in) :: l
80 : : complex(dp),intent(in) :: v1(l), v2(l)
81 : : ! *************************************************************************
82 : :
83 [ + + ]: 1288680 : complex_vector_product = sum(conjg(v1(:))*v2(:))
84 : :
85 : 9848 : end function complex_vector_product
86 : : !!***
87 : :
88 : : !!****f* m_gwls_utility/orthogonalize
89 : : !! NAME
90 : : !! orthogonalize
91 : : !!
92 : : !! FUNCTION
93 : : !! .
94 : : !!
95 : : !! INPUTS
96 : : !!
97 : : !! OUTPUT
98 : : !!
99 : : !! SOURCE
100 : :
101 : 1904 : subroutine orthogonalize(mpi_communicator, Hsize,Qsize,Xsize,Q,X)
102 : : !--------------------------------------------------------------------------
103 : : ! This function wraps the relevant LAPACK routines to perform
104 : : !
105 : : ! X = (1 - Q.Q^H ) . X
106 : : !
107 : : ! This essentially projects out the subspace in Q from X.
108 : : ! The array dimensions are:
109 : : !
110 : : ! Q [ Hsize, Qsize]
111 : : ! X [ Hsize, Xsize]
112 : : ! Y [ Hsize, Xsize]
113 : : !
114 : : ! Hsize means "dimension of the Hilbert space", so typically the number
115 : : ! of plane waves...
116 : : !--------------------------------------------------------------------------
117 : : integer, intent(in) :: mpi_communicator
118 : : integer, intent(in) :: Hsize, Qsize, Xsize
119 : : complex(dp),intent(in) :: Q(Hsize,Qsize)
120 : : complex(dp),intent(inout) :: X(Hsize,Xsize)
121 : 1904 : complex(dp),allocatable :: C(:,:)
122 : :
123 : : integer :: ierr
124 : :
125 : : ! *************************************************************************
126 : :
127 : :
128 [ - + + - : 7616 : ABI_MALLOC(C,(Qsize,Xsize))
- + + - +
- + - - +
- + ]
129 : :
130 : : ! Compute Q^dagger . X
131 : : call ZGEMM( 'C', & ! Hermitian conjugate the first array
132 : : 'N', & ! Leave second array as is
133 : : Qsize, & ! the number of rows of the matrix op( A )
134 : : Xsize, & ! the number of columns of the matrix op( B )
135 : : Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
136 : : cmplx_1, & ! alpha constant
137 : : Q, & ! matrix A
138 : : Hsize, & ! LDA
139 : : X, & ! matrix B
140 : : Hsize, & ! LDB
141 : : cmplx_0, & ! beta constant
142 : : C, & ! matrix C
143 : 1904 : Qsize) ! LDC
144 : :
145 : :
146 : 1904 : call xmpi_sum(C,mpi_communicator,ierr) ! sum on all processors working on FFT!
147 : :
148 : : ! Compute X - Q.(Q^dagger . X)
149 : : call ZGEMM( 'N', & ! Leave first array as is
150 : : 'N', & ! Leave second array as is
151 : : Hsize, & ! the number of rows of the matrix op( A )
152 : : Xsize, & ! the number of columns of the matrix op( B )
153 : : Qsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
154 : : -cmplx_1, & ! alpha constant
155 : : Q, & ! matrix A
156 : : Hsize, & ! LDA
157 : : C, & ! matrix B
158 : : Qsize, & ! LDB
159 : : cmplx_1, & ! beta constant
160 : : X, & ! matrix C
161 : 1904 : Hsize) ! LDC
162 : :
163 [ - + ]: 1904 : ABI_FREE(C)
164 : :
165 : 1904 : end subroutine orthogonalize
166 : : !!***
167 : :
168 : : !!****f* m_gwls_utility/driver_invert_positive_definite_hermitian_matrix
169 : : !! NAME
170 : : !! driver_invert_positive_definite_hermitian_matrix
171 : : !!
172 : : !! FUNCTION
173 : : !! .
174 : : !!
175 : : !! INPUTS
176 : : !!
177 : : !! OUTPUT
178 : : !!
179 : : !! SOURCE
180 : :
181 : 175 : subroutine driver_invert_positive_definite_hermitian_matrix(matrix,ldim)
182 : : !----------------------------------------------------------------------------------------------------
183 : : ! This is a utility-type subroutine, which encapsulates the many Lapack steps necessary
184 : : ! to invert a positive definite hermitian matrix, and returns the full inverse to avoid
185 : : ! errors!
186 : : !
187 : : ! The subroutine overwrites the input.
188 : : !----------------------------------------------------------------------------------------------------
189 : : integer , intent(in) :: ldim
190 : : complex(dp), intent(inout) :: matrix(ldim,ldim)
191 : :
192 : : ! local variables
193 : : integer :: i, j
194 : :
195 : :
196 : : integer :: info
197 : :
198 : : ! *************************************************************************
199 : :
200 : :
201 : :
202 : : ! First, peform a decomposition
203 : 175 : call abi_xpotrf( 'U', ldim, matrix, ldim, info )
204 : :
205 : : ! Second, inverse the matrix in the new format
206 : 175 : call zpotri( 'U', ldim,matrix, ldim, info )
207 : :
208 : :
209 : : ! Finally, properly symmetrise the matrix so that it is hermitian!
210 : : ! The upper triangular part of the matrix is correct; the lower triangular must be built
211 [ + + ]: 1571 : do j=1, ldim
212 [ + + ]: 6449 : do i=j+1, ldim
213 : 6274 : matrix(i,j)= conjg(matrix(j,i))
214 : : end do
215 : : end do
216 : :
217 : :
218 : 175 : end subroutine driver_invert_positive_definite_hermitian_matrix
219 : : !!***
220 : :
221 : : !!****f* m_gwls_utility/ritz_analysis_general
222 : : !! NAME
223 : : !! ritz_analysis_general
224 : : !!
225 : : !! FUNCTION
226 : : !! .
227 : : !!
228 : : !! INPUTS
229 : : !!
230 : : !! OUTPUT
231 : : !!
232 : : !! SOURCE
233 : :
234 : 2 : subroutine ritz_analysis_general(mpi_communicator,matrix_function,lmax,Hsize,Lbasis,eigenvalues)
235 : : !----------------------------------------------------------------------
236 : : !
237 : : ! This subroutine is mostly for testing purposes.
238 : : !
239 : : ! Given a matrix A_{NxN} which undergoes Lanczos analysis to yield
240 : : ! a trigonal matrix T_{k x k}, for k lanczos steps, it is useful to
241 : : ! test how well the eivenalues of T reproduce the eigenvalues of A.
242 : : !
243 : : ! Following Matrix computations by Golub and Van Loan, define
244 : : !
245 : : ! Q = [ | | ... | ], T = S^H.D.S, D diagonal, S unitary
246 : : ! [ q1 q2 ... qk]
247 : : ! [ | | ... | ]
248 : : !
249 : : ! Y = Q.S
250 : : !
251 : : ! If tridiagonalisation was taken all the way to k = N, we would expect
252 : : ! Y to contain the eigenvectors of A. It is useful to ask if, for a finite
253 : : ! k, Y already contains good approximations to eigenvectors, by computing
254 : : ! the Ritz residual,
255 : : !
256 : : ! Ri = || A.yi - di yi ||.
257 : : !
258 : : ! INPUTS:
259 : : ! matrix_function : the function which yields the action
260 : : ! of the implicit matrix on a vector
261 : : ! lmax : the total number of Lanczos steps
262 : : ! Hsize : the dimension of the Hilbert space
263 : : ! Lbasis : the Y matrix
264 : : ! eigenvalues : the computed approximate eigenvalues
265 : : ! of the matrix
266 : : !----------------------------------------------------------------------
267 : : interface
268 : : subroutine matrix_function(v_out,v_in,l)
269 : : use defs_basis
270 : :
271 : : integer, intent(in) :: l
272 : : complex(dp), intent(out) :: v_out(l)
273 : : complex(dp), intent(in) :: v_in(l)
274 : :
275 : : end subroutine matrix_function
276 : : end interface
277 : :
278 : :
279 : : integer, intent(in) :: Hsize, lmax , mpi_communicator
280 : : complex(dp), intent(in) :: Lbasis(Hsize,lmax)
281 : : real(dp), intent(in) :: eigenvalues(lmax)
282 : :
283 : :
284 : : ! local variables
285 : 2 : complex(dp),allocatable :: check_matrix(:,:)
286 : 2 : complex(dp),allocatable :: yl(:), rl(:), Ayl(:)
287 : :
288 : : real(dp) :: lambda_l
289 : : real(dp) :: check_norm
290 : : integer :: l, i
291 : :
292 : : real(dp) :: norm_ritz, norm_ritz_squared
293 : :
294 : : character(128) :: filename
295 : : logical :: file_exists
296 : : integer :: io_unit
297 : : integer :: ierr
298 : : integer :: mpi_rank
299 : :
300 : : logical :: head_node
301 : :
302 : : ! *************************************************************************
303 : :
304 [ - + + - : 8 : ABI_MALLOC(check_matrix,(lmax,lmax))
- + + - +
- + - - +
- + ]
305 [ - + + - : 6 : ABI_MALLOC(yl,(Hsize))
+ - + - -
+ - + ]
306 [ - + + - : 4 : ABI_MALLOC(rl,(Hsize))
- + - + ]
307 [ - + + - : 4 : ABI_MALLOC(Ayl,(Hsize))
- + - + ]
308 : :
309 : 2 : mpi_rank = xmpi_comm_rank(mpi_communicator)
310 : :
311 : 2 : head_node = mpi_rank == 0
312 : :
313 [ + - ]: 2 : if (head_node) then
314 : 2 : io_unit = get_unit()
315 : :
316 : 2 : i = 0
317 : 2 : file_exists = .true.
318 [ + + ]: 5 : do while (file_exists)
319 : 3 : i = i+1
320 : 3 : write(filename,'(A,I0.4,A)') "General_Ritz_Analisis_",i,".log"
321 : 5 : inquire(file=filename,exist=file_exists)
322 : : end do
323 : :
324 : :
325 : 2 : open(io_unit,file=filename,status=files_status_new)
326 : :
327 : 2 : write(io_unit,10) ''
328 : 2 : write(io_unit,10) '#===================================================================================================='
329 : 2 : write(io_unit,10) '# Entering ritz_analisis'
330 : 2 : write(io_unit,10) '# '
331 : 2 : write(io_unit,10) '# parameters'
332 : 2 : write(io_unit,16) '# Dimension of Hilbert space : ',Hsize
333 : 2 : write(io_unit,16) '# total number of Lanczos steps : ',lmax
334 : 2 : write(io_unit,10) '# '
335 : 2 : write(io_unit,10) '#===================================================================================================='
336 : 2 : write(io_unit,10) ''
337 : 2 : flush(io_unit)
338 : : end if
339 : :
340 : : ! NEVER use MATMUL! It stores temp arrays on stack, which kills executables compiled with intel!
341 : : ! check_matrix(:,:) = matmul(transpose(conjg(Lbasis)),Lbasis)
342 : : call ZGEMM( 'C', & ! Hermitian conjugate the first array
343 : : 'N', & ! Leave second array as is
344 : : lmax, & ! the number of rows of the matrix op( A )
345 : : lmax, & ! the number of columns of the matrix op( B )
346 : : Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
347 : : cmplx_1, & ! alpha constant
348 : : Lbasis, & ! matrix A
349 : : Hsize, & ! LDA
350 : : Lbasis, & ! matrix B
351 : : Hsize, & ! LDB
352 : : cmplx_0, & ! beta constant
353 : : check_matrix, & ! matrix C
354 : 2 : lmax) ! LDC
355 : :
356 : :
357 : 2 : call xmpi_sum(check_matrix,mpi_communicator,ierr) ! sum on all processors working on FFT!
358 : :
359 [ + + ]: 18 : do l = 1, lmax
360 : 18 : check_matrix(l,l) = check_matrix(l,l) - cmplx_1
361 : : end do
362 [ + + + + ]: 146 : check_norm = sqrt(sum(abs(check_matrix(:,:))**2))
363 : :
364 : :
365 [ + - ]: 2 : if (head_node) then
366 : 2 : write(io_unit,10) '#'
367 : 2 : write(io_unit,11) '# Is the basis orthonormal? || I - Y^H . Y || = ',check_norm
368 : 2 : write(io_unit,10) '#'
369 : 2 : flush(io_unit)
370 : :
371 : :
372 : : ! Compute Ritz norms
373 : 2 : write(io_unit,10) ''
374 : 2 : write(io_unit,10) '# Ritz analysis Lanczos Basis'
375 : 2 : write(io_unit,10) "# l lambda_l || R_l || "
376 : 2 : write(io_unit,10) '#===================================================================================================='
377 : 2 : flush(io_unit)
378 : : end if
379 : :
380 [ + + ]: 18 : do l = 1, lmax
381 : :
382 : 16 : lambda_l = eigenvalues(l)
383 : :
384 [ - + + + ]: 4144 : yl = Lbasis(:,l)
385 : 16 : call matrix_function(Ayl,yl,Hsize)
386 [ - + + + ]: 4144 : rl = Ayl - lambda_l*yl
387 : :
388 [ + + ]: 4128 : norm_ritz_squared = sum(abs(rl(:))**2)
389 : 16 : call xmpi_sum(norm_ritz_squared ,mpi_communicator,ierr)
390 : :
391 : 16 : norm_ritz = sqrt(norm_ritz_squared )
392 : :
393 : :
394 [ + - ]: 18 : if (head_node) write(io_unit,14) l, lambda_l, norm_ritz
395 : :
396 : : end do
397 : :
398 [ + - ]: 2 : if (head_node) then
399 : 2 : flush(io_unit)
400 : 2 : close(io_unit)
401 : : end if
402 : :
403 : :
404 [ - + ]: 2 : ABI_FREE(check_matrix)
405 [ - + ]: 2 : ABI_FREE(yl)
406 [ - + ]: 2 : ABI_FREE(rl)
407 [ - + ]: 2 : ABI_FREE(Ayl)
408 : :
409 : :
410 : : 10 format(A)
411 : : 11 format(A,ES24.16)
412 : : 14 format(I5,2(5X,F24.12))
413 : : 16 format(A,I5)
414 : :
415 : :
416 : 2 : end subroutine ritz_analysis_general
417 : : !!***
418 : :
419 : : end module m_gwls_utility
420 : :
421 : : !!***
|
