TABLE OF CONTENTS


ABINIT/m_gwls_utility [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_utility

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 
28 !---------------------------------------------------------------------
29 !  Utility modules, which are not directly related to the
30 !  problem to be solved.
31 !---------------------------------------------------------------------
32 
33 module m_gwls_utility
34 !----------------------------------------------------------------------------------------------------
35 ! This module contains useful functions and subroutines, which are not necessarily numerical
36 ! in nature, for timing, opening files, etc...
37 !----------------------------------------------------------------------------------------------------
38 
39 ! abinit modules
40 use defs_basis
41 use defs_datatypes
42 use m_abicore
43 use m_xmpi
44 
45 use m_io_tools, only : get_unit
46 
47 implicit none
48 private
49 
50 complex(dpc), public, parameter :: cmplx_i = (0.0_dp,1.0_dp)
51 complex(dpc), public, parameter :: cmplx_1 = (1.0_dp,0.0_dp)
52 complex(dpc), public, parameter :: cmplx_0 = (0.0_dp,0.0_dp)
53 
54 logical, public  :: master_debug
55 character(len=100), public :: files_status_new='new'
56 character(len=100), public :: files_status_old='unknown'
57 
58 public :: driver_invert_positive_definite_hermitian_matrix
59 public :: ritz_analysis_general, orthogonalize
60 public :: complex_vector_product

m_gwls_utility/complex_vector_product [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  complex_vector_product

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

 83 complex(dpc) function complex_vector_product(v1,v2,l)
 84 !--------------------------------------------------------------------------
 85 ! This function computes the vector product of two complex vectors.
 86 !--------------------------------------------------------------------------
 87 
 88 !This section has been created automatically by the script Abilint (TD).
 89 !Do not modify the following lines by hand.
 90 #undef ABI_FUNC
 91 #define ABI_FUNC 'complex_vector_product'
 92 !End of the abilint section
 93 
 94 implicit none
 95 
 96 integer,     intent(in)  :: l
 97 complex(dpc),intent(in)  :: v1(l), v2(l)
 98 
 99 ! *************************************************************************
100 
101 complex_vector_product = sum(conjg(v1(:))*v2(:))
102 
103 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

PARENTS

      gwls_DielectricArray,gwls_GenerateEpsilon

CHILDREN

      matrix_function,xmpi_sum,zgemm

SOURCE

222 subroutine driver_invert_positive_definite_hermitian_matrix(matrix,ldim)
223 !----------------------------------------------------------------------------------------------------
224 !        This is a utility-type subroutine, which encapsulates the many Lapack steps necessary
225 !        to invert a positive definite hermitian matrix, and returns the full inverse to avoid
226 !        errors!
227 !
228 !        The subroutine overwrites the input.
229 !----------------------------------------------------------------------------------------------------
230 
231 !This section has been created automatically by the script Abilint (TD).
232 !Do not modify the following lines by hand.
233 #undef ABI_FUNC
234 #define ABI_FUNC 'driver_invert_positive_definite_hermitian_matrix'
235 !End of the abilint section
236 
237 implicit none
238 
239 integer     , intent(in)    :: ldim
240 complex(dpc), intent(inout) :: matrix(ldim,ldim)
241 
242 ! local variables
243 integer      :: i, j
244 
245 
246 integer      :: info
247 
248 ! *************************************************************************
249 
250 
251 
252 ! First, peform a decomposition
253 call zpotrf( 'U', ldim,matrix, ldim, info )
254 
255 ! Second, inverse the matrix in the new format
256 call zpotri( 'U', ldim,matrix, ldim, info )
257 
258 
259 ! Finally,  properly symmetrise the matrix so that it is hermitian!
260 ! The upper triangular part of the matrix is correct; the lower triangular must be built
261 do j=1, ldim
262 do i=j+1, ldim
263 matrix(i,j)= conjg(matrix(j,i))
264 end do
265 end do
266 
267 
268 end subroutine driver_invert_positive_definite_hermitian_matrix

m_gwls_utility/orthogonalize [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  orthogonalize

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_GWlanczos

CHILDREN

      matrix_function,xmpi_sum,zgemm

SOURCE

125 subroutine orthogonalize(mpi_communicator, Hsize,Qsize,Xsize,Q,X)
126 !--------------------------------------------------------------------------
127 ! This function wraps the relevant LAPACK routines to perform
128 !
129 !                X = (1 - Q.Q^H ) . X
130 !
131 ! This essentially projects out  the subspace in Q from X.
132 ! The array dimensions are:
133 !
134 !                                Q [ Hsize, Qsize]
135 !                                X [ Hsize, Xsize]
136 !                                Y [ Hsize, Xsize]
137 !
138 !  Hsize means "dimension of the Hilbert space", so typically the number
139 !  of plane waves...
140 !--------------------------------------------------------------------------
141 
142 !This section has been created automatically by the script Abilint (TD).
143 !Do not modify the following lines by hand.
144 #undef ABI_FUNC
145 #define ABI_FUNC 'orthogonalize'
146 !End of the abilint section
147 
148 implicit none
149 
150 integer,     intent(in)  :: mpi_communicator
151 integer,     intent(in)  :: Hsize, Qsize, Xsize
152 complex(dpc),intent(in)  :: Q(Hsize,Qsize)
153 
154 complex(dpc),intent(inout)  :: X(Hsize,Xsize)
155 
156 complex(dpc),allocatable :: C(:,:)
157 
158 integer :: ierr
159 
160 ! *************************************************************************
161 
162 
163 ABI_ALLOCATE(C,(Qsize,Xsize))
164 
165 ! Compute Q^dagger . X
166 call ZGEMM(            'C',   & ! Hermitian conjugate the first array
167 'N',   & ! Leave second array as is
168 Qsize,   & ! the number of rows of the  matrix op( A )
169 Xsize,   & ! the number of columns of the  matrix op( B )
170 Hsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
171 cmplx_1,   & ! alpha constant
172 Q,   & ! matrix A
173 Hsize,   & ! LDA
174 X,   & ! matrix B
175 Hsize,   & ! LDB
176 cmplx_0,   & ! beta constant
177 C,   & ! matrix C
178 Qsize)     ! LDC
179 
180 
181 call xmpi_sum(C,mpi_communicator,ierr) ! sum on all processors working on FFT!
182 
183 ! Compute X - Q.(Q^dagger . X)
184 call ZGEMM(            'N',   & ! Leave first array as is
185 'N',   & ! Leave second array as is
186 Hsize,   & ! the number of rows of the  matrix op( A )
187 Xsize,   & ! the number of columns of the  matrix op( B )
188 Qsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
189 -cmplx_1,   & ! alpha constant
190 Q,   & ! matrix A
191 Hsize,   & ! LDA
192 C,   & ! matrix B
193 Qsize,   & ! LDB
194 cmplx_1,   & ! beta constant
195 X,   & ! matrix C
196 Hsize)     ! LDC
197 
198 ABI_DEALLOCATE(C)
199 
200 end subroutine orthogonalize

m_gwls_utility/ritz_analysis_general [ Functions ]

[ Top ] [ m_gwls_utility ] [ Functions ]

NAME

  ritz_analysis_general

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputePoles,gwls_GenerateEpsilon

CHILDREN

      matrix_function,xmpi_sum,zgemm

SOURCE

290 subroutine ritz_analysis_general(mpi_communicator,matrix_function,lmax,Hsize,Lbasis,eigenvalues)
291 !----------------------------------------------------------------------
292 !
293 ! This subroutine is mostly for testing purposes.
294 !
295 ! Given a matrix A_{NxN} which undergoes Lanczos analysis to yield
296 ! a trigonal matrix T_{k x k}, for k lanczos steps, it is useful to
297 ! test how well the eivenalues of T reproduce the eigenvalues of A.
298 !
299 ! Following Matrix computations by Golub and Van Loan, define
300 !
301 !                Q = [ |  |  ...  | ],        T = S^H.D.S,  D diagonal, S unitary
302 !                    [ q1 q2 ...  qk]
303 !                    [ |  |  ...  | ]
304 !
305 !                Y = Q.S
306 !
307 ! If tridiagonalisation was taken all the way to k = N, we would expect
308 ! Y to contain the eigenvectors of A. It is useful to ask if, for a finite
309 ! k, Y already contains good approximations to eigenvectors, by computing
310 ! the Ritz residual,
311 !
312 !        Ri = || A.yi - di yi ||.
313 !
314 !        INPUTS:
315 !                     matrix_function    : the function which yields the action
316 !                                        of the implicit matrix on a vector
317 !                        lmax            : the total number of Lanczos steps
318 !                        Hsize           : the dimension of the Hilbert space
319 !                        Lbasis          : the Y matrix
320 !                       eigenvalues      : the computed approximate eigenvalues
321 !                                        of the matrix
322 !----------------------------------------------------------------------
323 
324 !This section has been created automatically by the script Abilint (TD).
325 !Do not modify the following lines by hand.
326 #undef ABI_FUNC
327 #define ABI_FUNC 'ritz_analysis_general'
328 !End of the abilint section
329 
330 implicit none
331 
332 interface
333   subroutine matrix_function(v_out,v_in,l)
334   use defs_basis
335 
336   integer,     intent(in)  :: l
337   complex(dp), intent(out) :: v_out(l)
338   complex(dp), intent(in)  :: v_in(l)
339 
340   end subroutine matrix_function
341 end interface
342 
343 
344 integer,      intent(in)    :: Hsize, lmax , mpi_communicator
345 complex(dpc), intent(in)    :: Lbasis(Hsize,lmax)
346 real(dp),     intent(in)    :: eigenvalues(lmax)
347 
348 
349 ! local variables
350 complex(dpc),allocatable :: check_matrix(:,:)
351 complex(dpc),allocatable :: yl(:), rl(:), Ayl(:)
352 
353 real(dp)     :: lambda_l
354 real(dp)     :: check_norm
355 integer      :: l, i
356 
357 real(dp)     :: norm_ritz, norm_ritz_squared
358 
359 character(128) :: filename
360 logical        :: file_exists
361 integer        :: io_unit
362 integer        :: ierr
363 integer        :: mpi_rank
364 
365 logical        :: head_node
366 
367 ! *************************************************************************
368 
369 ABI_ALLOCATE(check_matrix,(lmax,lmax))
370 ABI_ALLOCATE(yl,(Hsize))
371 ABI_ALLOCATE(rl,(Hsize))
372 ABI_ALLOCATE(Ayl,(Hsize))
373 
374 mpi_rank  = xmpi_comm_rank(mpi_communicator)
375 
376 head_node = mpi_rank == 0
377 
378 if (head_node) then
379   io_unit = get_unit()
380 
381   i = 0
382   file_exists = .true.
383   do while (file_exists)
384   i = i+1
385   write(filename,'(A,I0.4,A)') "General_Ritz_Analisis_",i,".log"
386   inquire(file=filename,exist=file_exists)
387   end do
388 
389 
390   open(io_unit,file=filename,status=files_status_new)
391 
392   write(io_unit,10) ''
393   write(io_unit,10) '#===================================================================================================='
394   write(io_unit,10) '# Entering ritz_analisis'
395   write(io_unit,10) '# '
396   write(io_unit,10) '#  parameters'
397   write(io_unit,16) '#          Dimension of Hilbert space    : ',Hsize
398   write(io_unit,16) '#          total number of Lanczos steps : ',lmax
399   write(io_unit,10) '# '
400   write(io_unit,10) '#===================================================================================================='
401   write(io_unit,10) ''
402   flush(io_unit)
403 end if
404 
405 ! NEVER use MATMUL! It stores temp arrays on stack, which kills executables compiled with intel!
406 ! check_matrix(:,:) = matmul(transpose(conjg(Lbasis)),Lbasis)
407 call ZGEMM(            'C',   & ! Hermitian conjugate the first array
408 'N',   & ! Leave second array as is
409 lmax,   & ! the number of rows of the  matrix op( A )
410 lmax,   & ! the number of columns of the  matrix op( B )
411 Hsize,   & ! the number of columns of the  matrix op( A ) == rows of matrix op( B )
412 cmplx_1,   & ! alpha constant
413 Lbasis,   & ! matrix A
414 Hsize,   & ! LDA
415 Lbasis,   & ! matrix B
416 Hsize,   & ! LDB
417 cmplx_0,   & ! beta constant
418 check_matrix,   & ! matrix C
419 lmax)     ! LDC
420 
421 
422 call xmpi_sum(check_matrix,mpi_communicator,ierr) ! sum on all processors working on FFT!
423 
424 do l = 1, lmax
425 check_matrix(l,l) =  check_matrix(l,l) - cmplx_1
426 end do
427 check_norm = sqrt(sum(abs(check_matrix(:,:))**2))
428 
429 
430 if (head_node) then
431   write(io_unit,10) '#'
432   write(io_unit,11) '#  Is the basis orthonormal?   || I - Y^H . Y || = ',check_norm
433   write(io_unit,10) '#'
434   flush(io_unit)
435 
436 
437   ! Compute Ritz norms
438   write(io_unit,10) ''
439   write(io_unit,10) '#  Ritz analysis Lanczos Basis'
440   write(io_unit,10) "#   l               lambda_l                       || R_l ||                                         "
441   write(io_unit,10) '#===================================================================================================='
442   flush(io_unit)
443 end if
444 
445 do l = 1, lmax
446 
447 lambda_l = eigenvalues(l)
448 
449 yl       = Lbasis(:,l)
450 call matrix_function(Ayl,yl,Hsize)
451 rl       = Ayl - lambda_l*yl
452 
453 norm_ritz_squared = sum(abs(rl(:))**2)
454 call xmpi_sum(norm_ritz_squared ,mpi_communicator,ierr)
455 
456 norm_ritz = sqrt(norm_ritz_squared )
457 
458 
459 if (head_node) write(io_unit,14) l, lambda_l, norm_ritz
460 
461 end do
462 
463 if (head_node) then
464   flush(io_unit)
465   close(io_unit)
466 end if
467 
468 
469 ABI_DEALLOCATE(check_matrix)
470 ABI_DEALLOCATE(yl)
471 ABI_DEALLOCATE(rl)
472 ABI_DEALLOCATE(Ayl)
473 
474 
475 10 format(A)
476 11 format(A,ES24.16)
477 14 format(I5,2(5X,F24.12))
478 16 format(A,I5)
479 
480 
481 end subroutine ritz_analysis_general