TABLE OF CONTENTS


ABINIT/m_gwls_utility [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_utility

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2022 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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_io_tools, only : get_unit
40 
41 implicit none
42 
43 private
44 
45 complex(dpc), public, parameter :: cmplx_i = (0.0_dp,1.0_dp)
46 complex(dpc), public, parameter :: cmplx_1 = (1.0_dp,0.0_dp)
47 complex(dpc), public, parameter :: cmplx_0 = (0.0_dp,0.0_dp)
48 
49 logical, public  :: master_debug
50 character(len=100), public :: files_status_new='new'
51 character(len=100), public :: files_status_old='unknown'
52 
53 public :: driver_invert_positive_definite_hermitian_matrix
54 public :: ritz_analysis_general, orthogonalize
55 public :: complex_vector_product

m_gwls_utility/complex_vector_product [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  complex_vector_product

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

74 complex(dpc) function complex_vector_product(v1,v2,l)
75 !--------------------------------------------------------------------------
76 ! This function computes the vector product of two complex vectors.
77 !--------------------------------------------------------------------------
78 implicit none
79 
80 integer,     intent(in)  :: l
81 complex(dpc),intent(in)  :: v1(l), v2(l)
82 
83 ! *************************************************************************
84 
85 complex_vector_product = sum(conjg(v1(:))*v2(:))
86 
87 end function complex_vector_product

m_gwls_utility/driver_invert_positive_definite_hermitian_matrix [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  driver_invert_positive_definite_hermitian_matrix

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

187 subroutine driver_invert_positive_definite_hermitian_matrix(matrix,ldim)
188 !----------------------------------------------------------------------------------------------------
189 !        This is a utility-type subroutine, which encapsulates the many Lapack steps necessary
190 !        to invert a positive definite hermitian matrix, and returns the full inverse to avoid
191 !        errors!
192 !
193 !        The subroutine overwrites the input.
194 !----------------------------------------------------------------------------------------------------
195 implicit none
196 
197 integer     , intent(in)    :: ldim
198 complex(dpc), intent(inout) :: matrix(ldim,ldim)
199 
200 ! local variables
201 integer      :: i, j
202 
203 
204 integer      :: info
205 
206 ! *************************************************************************
207 
208 
209 
210 ! First, peform a decomposition
211 call zpotrf( 'U', ldim,matrix, ldim, info )
212 
213 ! Second, inverse the matrix in the new format
214 call zpotri( 'U', ldim,matrix, ldim, info )
215 
216 
217 ! Finally,  properly symmetrise the matrix so that it is hermitian!
218 ! The upper triangular part of the matrix is correct; the lower triangular must be built
219 do j=1, ldim
220 do i=j+1, ldim
221 matrix(i,j)= conjg(matrix(j,i))
222 end do
223 end do
224 
225 
226 end subroutine driver_invert_positive_definite_hermitian_matrix

m_gwls_utility/orthogonalize [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  orthogonalize

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

103 subroutine orthogonalize(mpi_communicator, Hsize,Qsize,Xsize,Q,X)
104 !--------------------------------------------------------------------------
105 ! This function wraps the relevant LAPACK routines to perform
106 !
107 !                X = (1 - Q.Q^H ) . X
108 !
109 ! This essentially projects out  the subspace in Q from X.
110 ! The array dimensions are:
111 !
112 !                                Q [ Hsize, Qsize]
113 !                                X [ Hsize, Xsize]
114 !                                Y [ Hsize, Xsize]
115 !
116 !  Hsize means "dimension of the Hilbert space", so typically the number
117 !  of plane waves...
118 !--------------------------------------------------------------------------
119 implicit none
120 
121 integer,     intent(in)  :: mpi_communicator
122 integer,     intent(in)  :: Hsize, Qsize, Xsize
123 complex(dpc),intent(in)  :: Q(Hsize,Qsize)
124 
125 complex(dpc),intent(inout)  :: X(Hsize,Xsize)
126 
127 complex(dpc),allocatable :: C(:,:)
128 
129 integer :: ierr
130 
131 ! *************************************************************************
132 
133 
134 ABI_MALLOC(C,(Qsize,Xsize))
135 
136 ! Compute Q^dagger . X
137 call ZGEMM(            'C',   & ! Hermitian conjugate the first array
138 'N',   & ! Leave second array as is
139 Qsize,   & ! the number of rows of the  matrix op( A )
140 Xsize,   & ! the number of columns of the  matrix op( B )
141 Hsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
142 cmplx_1,   & ! alpha constant
143 Q,   & ! matrix A
144 Hsize,   & ! LDA
145 X,   & ! matrix B
146 Hsize,   & ! LDB
147 cmplx_0,   & ! beta constant
148 C,   & ! matrix C
149 Qsize)     ! LDC
150 
151 
152 call xmpi_sum(C,mpi_communicator,ierr) ! sum on all processors working on FFT!
153 
154 ! Compute X - Q.(Q^dagger . X)
155 call ZGEMM(            'N',   & ! Leave first array as is
156 'N',   & ! Leave second array as is
157 Hsize,   & ! the number of rows of the  matrix op( A )
158 Xsize,   & ! the number of columns of the  matrix op( B )
159 Qsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
160 -cmplx_1,   & ! alpha constant
161 Q,   & ! matrix A
162 Hsize,   & ! LDA
163 C,   & ! matrix B
164 Qsize,   & ! LDB
165 cmplx_1,   & ! beta constant
166 X,   & ! matrix C
167 Hsize)     ! LDC
168 
169 ABI_FREE(C)
170 
171 end subroutine orthogonalize

m_gwls_utility/ritz_analysis_general [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  ritz_analysis_general

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

242 subroutine ritz_analysis_general(mpi_communicator,matrix_function,lmax,Hsize,Lbasis,eigenvalues)
243 !----------------------------------------------------------------------
244 !
245 ! This subroutine is mostly for testing purposes.
246 !
247 ! Given a matrix A_{NxN} which undergoes Lanczos analysis to yield
248 ! a trigonal matrix T_{k x k}, for k lanczos steps, it is useful to
249 ! test how well the eivenalues of T reproduce the eigenvalues of A.
250 !
251 ! Following Matrix computations by Golub and Van Loan, define
252 !
253 !                Q = [ |  |  ...  | ],        T = S^H.D.S,  D diagonal, S unitary
254 !                    [ q1 q2 ...  qk]
255 !                    [ |  |  ...  | ]
256 !
257 !                Y = Q.S
258 !
259 ! If tridiagonalisation was taken all the way to k = N, we would expect
260 ! Y to contain the eigenvectors of A. It is useful to ask if, for a finite
261 ! k, Y already contains good approximations to eigenvectors, by computing
262 ! the Ritz residual,
263 !
264 !        Ri = || A.yi - di yi ||.
265 !
266 !        INPUTS:
267 !                     matrix_function    : the function which yields the action
268 !                                        of the implicit matrix on a vector
269 !                        lmax            : the total number of Lanczos steps
270 !                        Hsize           : the dimension of the Hilbert space
271 !                        Lbasis          : the Y matrix
272 !                       eigenvalues      : the computed approximate eigenvalues
273 !                                        of the matrix
274 !----------------------------------------------------------------------
275 implicit none
276 
277 interface
278   subroutine matrix_function(v_out,v_in,l)
279   use defs_basis
280 
281   integer,     intent(in)  :: l
282   complex(dp), intent(out) :: v_out(l)
283   complex(dp), intent(in)  :: v_in(l)
284 
285   end subroutine matrix_function
286 end interface
287 
288 
289 integer,      intent(in)    :: Hsize, lmax , mpi_communicator
290 complex(dpc), intent(in)    :: Lbasis(Hsize,lmax)
291 real(dp),     intent(in)    :: eigenvalues(lmax)
292 
293 
294 ! local variables
295 complex(dpc),allocatable :: check_matrix(:,:)
296 complex(dpc),allocatable :: yl(:), rl(:), Ayl(:)
297 
298 real(dp)     :: lambda_l
299 real(dp)     :: check_norm
300 integer      :: l, i
301 
302 real(dp)     :: norm_ritz, norm_ritz_squared
303 
304 character(128) :: filename
305 logical        :: file_exists
306 integer        :: io_unit
307 integer        :: ierr
308 integer        :: mpi_rank
309 
310 logical        :: head_node
311 
312 ! *************************************************************************
313 
314 ABI_MALLOC(check_matrix,(lmax,lmax))
315 ABI_MALLOC(yl,(Hsize))
316 ABI_MALLOC(rl,(Hsize))
317 ABI_MALLOC(Ayl,(Hsize))
318 
319 mpi_rank  = xmpi_comm_rank(mpi_communicator)
320 
321 head_node = mpi_rank == 0
322 
323 if (head_node) then
324   io_unit = get_unit()
325 
326   i = 0
327   file_exists = .true.
328   do while (file_exists)
329   i = i+1
330   write(filename,'(A,I0.4,A)') "General_Ritz_Analisis_",i,".log"
331   inquire(file=filename,exist=file_exists)
332   end do
333 
334 
335   open(io_unit,file=filename,status=files_status_new)
336 
337   write(io_unit,10) ''
338   write(io_unit,10) '#===================================================================================================='
339   write(io_unit,10) '# Entering ritz_analisis'
340   write(io_unit,10) '# '
341   write(io_unit,10) '#  parameters'
342   write(io_unit,16) '#          Dimension of Hilbert space    : ',Hsize
343   write(io_unit,16) '#          total number of Lanczos steps : ',lmax
344   write(io_unit,10) '# '
345   write(io_unit,10) '#===================================================================================================='
346   write(io_unit,10) ''
347   flush(io_unit)
348 end if
349 
350 ! NEVER use MATMUL! It stores temp arrays on stack, which kills executables compiled with intel!
351 ! check_matrix(:,:) = matmul(transpose(conjg(Lbasis)),Lbasis)
352 call ZGEMM(            'C',   & ! Hermitian conjugate the first array
353 'N',   & ! Leave second array as is
354 lmax,   & ! the number of rows of the  matrix op( A )
355 lmax,   & ! the number of columns of the  matrix op( B )
356 Hsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
357 cmplx_1,   & ! alpha constant
358 Lbasis,   & ! matrix A
359 Hsize,   & ! LDA
360 Lbasis,   & ! matrix B
361 Hsize,   & ! LDB
362 cmplx_0,   & ! beta constant
363 check_matrix,   & ! matrix C
364 lmax)     ! LDC
365 
366 
367 call xmpi_sum(check_matrix,mpi_communicator,ierr) ! sum on all processors working on FFT!
368 
369 do l = 1, lmax
370 check_matrix(l,l) =  check_matrix(l,l) - cmplx_1
371 end do
372 check_norm = sqrt(sum(abs(check_matrix(:,:))**2))
373 
374 
375 if (head_node) then
376   write(io_unit,10) '#'
377   write(io_unit,11) '#  Is the basis orthonormal?   || I - Y^H . Y || = ',check_norm
378   write(io_unit,10) '#'
379   flush(io_unit)
380 
381 
382   ! Compute Ritz norms
383   write(io_unit,10) ''
384   write(io_unit,10) '#  Ritz analysis Lanczos Basis'
385   write(io_unit,10) "#   l               lambda_l                       || R_l ||                                         "
386   write(io_unit,10) '#===================================================================================================='
387   flush(io_unit)
388 end if
389 
390 do l = 1, lmax
391 
392 lambda_l = eigenvalues(l)
393 
394 yl       = Lbasis(:,l)
395 call matrix_function(Ayl,yl,Hsize)
396 rl       = Ayl - lambda_l*yl
397 
398 norm_ritz_squared = sum(abs(rl(:))**2)
399 call xmpi_sum(norm_ritz_squared ,mpi_communicator,ierr)
400 
401 norm_ritz = sqrt(norm_ritz_squared )
402 
403 
404 if (head_node) write(io_unit,14) l, lambda_l, norm_ritz
405 
406 end do
407 
408 if (head_node) then
409   flush(io_unit)
410   close(io_unit)
411 end if
412 
413 
414 ABI_FREE(check_matrix)
415 ABI_FREE(yl)
416 ABI_FREE(rl)
417 ABI_FREE(Ayl)
418 
419 
420 10 format(A)
421 11 format(A,ES24.16)
422 14 format(I5,2(5X,F24.12))
423 16 format(A,I5)
424 
425 
426 end subroutine ritz_analysis_general