TABLE OF CONTENTS
- ABINIT/m_gwls_utility
- m_gwls_utility/complex_vector_product
- m_gwls_utility/driver_invert_positive_definite_hermitian_matrix
- m_gwls_utility/orthogonalize
- m_gwls_utility/ritz_analysis_general
ABINIT/m_gwls_utility [ 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