diago_lanczos.f90 2.96 KB
Newer Older
Marco Govoni's avatar
Marco Govoni committed
1
!
Marco Govoni's avatar
Marco Govoni committed
2
! Copyright (C) 2015-2021 M. Govoni 
Marco Govoni's avatar
Marco Govoni committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file: 
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE diago_lanczos( bnorm, diago, subdiago, braket, ndim ) 
  !-----------------------------------------------------------------------
  !
  USE kinds,               ONLY : DP
  USE westcom,             ONLY : n_lanczos
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  REAL(DP),INTENT(IN) :: bnorm
  REAL(DP),INTENT(INOUT) :: diago(n_lanczos)
  REAL(DP),INTENT(INOUT) :: subdiago(n_lanczos-1)
  REAL(DP),INTENT(INOUT) :: braket(ndim,n_lanczos)
  INTEGER,INTENT(IN) :: ndim
  !
  ! Workspace
  !
  INTEGER :: info,il,ip
  REAL(DP),ALLOCATABLE :: z(:,:)
  REAL(DP),ALLOCATABLE :: work(:)
  REAL(DP),ALLOCATABLE :: braket_tilde(:,:)
  !
  ALLOCATE( z(n_lanczos,n_lanczos), work(MAX(1,2*n_lanczos-2)), braket_tilde(ndim,n_lanczos) )
  !
  ! Diagonalize
  !
  CALL DSTEQR ('I', n_lanczos, diago, subdiago, z, n_lanczos, work, info)
  !
  ! braket_tilde(ip,il) = brak(ip,:) * z(:,il)
  !  
  CALL DGEMM('N','N',ndim,n_lanczos,n_lanczos,bnorm,braket,ndim,z,n_lanczos,0._DP,braket_tilde,ndim)
  ! 
  ! braket_tilde(ip,il) = braket_tilde(ip,il) * z(1,il)
  !
  DO il=1,n_lanczos
     DO ip = 1, ndim
        braket(ip,il) = braket_tilde(ip,il) * z(1,il)
     ENDDO
  ENDDO 
  !
  DEALLOCATE( z, work, braket_tilde )
  !
END SUBROUTINE
!
!-----------------------------------------------------------------------
SUBROUTINE diago_lanczos_complex( bnorm, diago, subdiago, braket, ndim ) 
  !-----------------------------------------------------------------------
  !
  USE kinds,               ONLY : DP
  USE westcom,             ONLY : n_lanczos
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  REAL(DP),INTENT(IN) :: bnorm
  REAL(DP),INTENT(INOUT) :: diago(n_lanczos)
  REAL(DP),INTENT(INOUT) :: subdiago(n_lanczos-1)
  COMPLEX(DP),INTENT(INOUT) :: braket(ndim,n_lanczos)
  INTEGER,INTENT(IN) :: ndim
  !
  ! Workspace
  !
  INTEGER :: info,il,ip
  COMPLEX(DP) :: znorm, zzero
  COMPLEX(DP),ALLOCATABLE :: z(:,:)
  REAL(DP),ALLOCATABLE :: work(:)
  COMPLEX(DP),ALLOCATABLE :: braket_tilde(:,:)
  !
  znorm = CMPLX( bnorm, 0._DP, KIND=DP )
  zzero = ( 0._DP, 0._DP ) 
  !
  ALLOCATE( z(n_lanczos,n_lanczos), work(MAX(1,2*n_lanczos-2)), braket_tilde(ndim,n_lanczos) )
  !
  ! Diagonalize
  !
  CALL ZSTEQR ('I', n_lanczos, diago, subdiago, z, n_lanczos, work, info)
  !
  ! braket_tilde(ip,il) = brak(ip,:) * z(:,il)
  !  
  CALL ZGEMM('N','N',ndim,n_lanczos,n_lanczos,znorm,braket,ndim,z,n_lanczos,zzero,braket_tilde,ndim)
  ! 
  ! braket_tilde(ip,il) = braket_tilde(ip,il) * z(1,il)
  !
  DO il=1,n_lanczos
     DO ip = 1, ndim
        braket(ip,il) = braket_tilde(ip,il) * z(1,il)
     ENDDO
  ENDDO 
  !
  DEALLOCATE( z, work, braket_tilde )
  !
END SUBROUTINE