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