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-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_profiling_abi 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