rdiaghg.f90 9.05 KB
Newer Older
1
!
2
! Copyright (C) 2003-2013 Quantum ESPRESSO group
3 4 5 6 7 8 9 10
! 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 .
!
!
!----------------------------------------------------------------------------
SUBROUTINE rdiaghg( n, m, h, s, ldh, e, v )
11
  !----------------------------------------------------------------------------
12
  ! ... Hv=eSv, with H symmetric matrix, S overlap matrix.
13
  ! ... On output both matrix are unchanged
14
  !
15
  ! ... LAPACK version - uses both DSYGV and DSYGVX
16
  !
17
  USE la_param,          ONLY : DP
18
  USE mp,                ONLY : mp_bcast
19
  USE mp_bands_util,     ONLY : me_bgrp, root_bgrp, intra_bgrp_comm
20
  !
21
  IMPLICIT NONE
22
  !
23
  INTEGER, INTENT(IN) :: n, m, ldh
24 25 26
    ! dimension of the matrix to be diagonalized
    ! number of eigenstates to be calculated
    ! leading dimension of h, as declared in the calling pgm unit
27
  REAL(DP), INTENT(INOUT) :: h(ldh,n), s(ldh,n)
28 29
    ! matrix to be diagonalized
    ! overlap matrix
30
  !
31
  REAL(DP), INTENT(OUT) :: e(n)
32
    ! eigenvalues
33
  REAL(DP), INTENT(OUT) :: v(ldh,m)
34
    ! eigenvectors (column-wise)
35
  !
36
  INTEGER               :: lwork, nb, mm, info, i, j
37
    ! mm = number of calculated eigenvectors
38
  REAL(DP)              :: abstol
39 40
  REAL(DP), PARAMETER   :: one = 1_DP
  REAL(DP), PARAMETER   :: zero = 0_DP
41
  INTEGER,  ALLOCATABLE :: iwork(:), ifail(:)
42
  REAL(DP), ALLOCATABLE :: work(:), sdiag(:), hdiag(:)
43
  LOGICAL               :: all_eigenvalues
44
  INTEGER,  EXTERNAL    :: ILAENV
45
    ! ILAENV returns optimal block size "nb"
46
  !
47
  CALL start_clock( 'rdiaghg' )
48
  !
49 50
  ! ... only the first processor diagonalize the matrix
  !
51
  IF ( me_bgrp == root_bgrp ) THEN
52
     !
53
     ! ... save the diagonal of input S (it will be overwritten)
54
     !
55 56 57 58
     ALLOCATE( sdiag( n ) )
     DO i = 1, n
        sdiag(i) = s(i,i)
     END DO
59
     !
60
     all_eigenvalues = ( m == n )
61
     !
62
     ! ... check for optimal block size
63
     !
64
     nb = ILAENV( 1, 'DSYTRD', 'U', n, -1, -1, -1 )
65
     !
66
     IF ( nb < 5 .OR. nb >= n ) THEN
67
        !
68
        lwork = 8*n
69
        !
70
     ELSE
71
        !
72
        lwork = ( nb + 3 )*n
73 74 75
        !
     END IF
     !
76
     ALLOCATE( work( lwork ) )
77
     !
78
     IF ( all_eigenvalues ) THEN
79
        !
80
        ! ... calculate all eigenvalues
81
        !
82 83 84 85 86
        !$omp parallel do
        do i =1, n
           v(1:ldh,i) = h(1:ldh,i)
        end do
        !$omp end parallel do
87
        !
88
#if defined (__ESSL)
89
        !
90
        ! ... there is a name conflict between essl and lapack ...
91
        !
92
        CALL DSYGV( 1, v, ldh, s, ldh, e, v, ldh, n, work, lwork )
93
        !
94 95 96 97
        info = 0
#else
        CALL DSYGV( 1, 'V', 'U', n, v, ldh, s, ldh, e, work, lwork, info )
#endif
98
        !
99
     ELSE
100
        !
101
        ! ... calculate only m lowest eigenvalues
102
        !
103 104
        ALLOCATE( iwork( 5*n ) )
        ALLOCATE( ifail( n ) )
105
        !
106
        ! ... save the diagonal of input H (it will be overwritten)
107
        !
108
        ALLOCATE( hdiag( n ) )
109
        DO i = 1, n
110
           hdiag(i) = h(i,i)
111
        END DO
112
        !
113 114
        abstol = 0.D0
       ! abstol = 2.D0*DLAMCH( 'S' )
115
        !
116 117 118
        CALL DSYGVX( 1, 'V', 'I', 'U', n, h, ldh, s, ldh, &
                     0.D0, 0.D0, 1, m, abstol, mm, e, v, ldh, &
                     work, lwork, iwork, ifail, info )
119
        !
120 121
        DEALLOCATE( ifail )
        DEALLOCATE( iwork )
122
        !
123
        ! ... restore input H matrix from saved diagonal and lower triangle
124
        !
125
        !$omp parallel do
126
        DO i = 1, n
127
           h(i,i) = hdiag(i)
128
           DO j = i + 1, n
129
              h(i,j) = h(j,i)
130 131
           END DO
           DO j = n + 1, ldh
132
              h(j,i) = 0.0_DP
133 134
           END DO
        END DO
135
        !$omp end parallel do
136
        !
137
        DEALLOCATE( hdiag )
138
        !
139
     END IF
140
     !
141 142
     DEALLOCATE( work )
     !
143 144 145 146 147 148 149
     IF ( info > n ) THEN
        CALL errore( 'rdiaghg', 'S matrix not positive definite', ABS( info ) )
     ELSE IF ( info > 0 ) THEN
        CALL errore( 'rdiaghg', 'eigenvectors failed to converge', ABS( info ) )
     ELSE IF ( info < 0 ) THEN
        CALL errore( 'rdiaghg', 'incorrect call to DSYGV*', ABS( info ) )
     END IF
150
     
151 152
     ! ... restore input S matrix from saved diagonal and lower triangle
     !
153
     !$omp parallel do
154 155 156 157 158 159 160 161 162
     DO i = 1, n
        s(i,i) = sdiag(i)
        DO j = i + 1, n
           s(i,j) = s(j,i)
        END DO
        DO j = n + 1, ldh
           s(j,i) = 0.0_DP
        END DO
     END DO
163
     !$omp end parallel do
164
     !
165
     DEALLOCATE( sdiag )
166
     !
167
  END IF
168
  !
169 170
  ! ... broadcast eigenvectors and eigenvalues to all other processors
  !
171 172
  CALL mp_bcast( e, root_bgrp, intra_bgrp_comm )
  CALL mp_bcast( v, root_bgrp, intra_bgrp_comm )
173
  !
174
  CALL stop_clock( 'rdiaghg' )
175 176 177
  !
  RETURN
  !
178
END SUBROUTINE rdiaghg
179 180 181 182 183 184 185 186 187
!
!----------------------------------------------------------------------------
SUBROUTINE prdiaghg( n, h, s, ldh, e, v, desc )
  !----------------------------------------------------------------------------
  !
  ! ... calculates eigenvalues and eigenvectors of the generalized problem
  ! ... Hv=eSv, with H symmetric matrix, S overlap matrix.
  ! ... On output both matrix are unchanged
  !
188
  ! ... Parallel version with full data distribution
189
  !
190
  USE la_param,          ONLY : DP
191 192
  USE mp,                ONLY : mp_bcast
  USE descriptors,       ONLY : la_descriptor
193
  USE mp_diag,           ONLY : ortho_parent_comm
194
#if defined __SCALAPACK
195 196
  USE mp_diag,           ONLY : ortho_cntx, me_blacs, np_ortho, me_ortho, ortho_comm
  USE dspev_module,      ONLY : pdsyevd_drv
197
#endif
198 199 200 201 202 203 204
  !
  !
  IMPLICIT NONE
  !
  INTEGER, INTENT(IN) :: n, ldh
    ! dimension of the matrix to be diagonalized and number of eigenstates to be calculated
    ! leading dimension of h, as declared in the calling pgm unit
cavazzon's avatar
cavazzon committed
205
  REAL(DP), INTENT(INOUT) :: h(ldh,ldh), s(ldh,ldh)
206 207 208 209 210
    ! matrix to be diagonalized
    ! overlap matrix
  !
  REAL(DP), INTENT(OUT) :: e(n)
    ! eigenvalues
cavazzon's avatar
cavazzon committed
211
  REAL(DP), INTENT(OUT) :: v(ldh,ldh)
212
    ! eigenvectors (column-wise)
213
  TYPE(la_descriptor), INTENT(IN) :: desc
214
  !
215
  INTEGER, PARAMETER    :: root = 0
216
  INTEGER               :: nx
217 218 219 220 221
    ! local block size
  REAL(DP), PARAMETER   :: one = 1_DP
  REAL(DP), PARAMETER   :: zero = 0_DP
  REAL(DP), ALLOCATABLE :: hh(:,:)
  REAL(DP), ALLOCATABLE :: ss(:,:)
222
#if defined(__SCALAPACK)
223 224
  INTEGER     :: desch( 16 ), info
#endif
225
  INTEGER               :: i
226
  !
227
  CALL start_clock( 'rdiaghg' )
228
  !
229
  IF( desc%active_node > 0 ) THEN
230
     !
231
     nx   = desc%nrcx
232 233 234 235 236 237 238
     !
     IF( nx /= ldh ) &
        CALL errore(" prdiaghg ", " inconsistent leading dimension ", ldh )
     !
     ALLOCATE( hh( nx, nx ) )
     ALLOCATE( ss( nx, nx ) )
     !
239 240 241 242 243 244
     !$omp parallel do
     do i=1,nx
        hh(1:nx,i) = h(1:nx,i)
        ss(1:nx,i) = s(1:nx,i)
     end do
     !$omp end parallel do
245 246 247
     !
  END IF
  !
248
  CALL start_clock( 'rdiaghg:choldc' )
249 250 251
  !
  ! ... Cholesky decomposition of s ( L is stored in s )
  !
252
  IF( desc%active_node > 0 ) THEN
253
     !
254
#if defined(__SCALAPACK)
255
     CALL descinit( desch, n, n, desc%nrcx, desc%nrcx, 0, 0, ortho_cntx, SIZE( hh, 1 ) , info )
256
  
257
     IF( info /= 0 ) CALL errore( ' rdiaghg ', ' descinit ', ABS( info ) )
258 259
#endif
     !
260
#if defined(__SCALAPACK)
261 262 263
     CALL PDPOTRF( 'L', n, ss, 1, 1, desch, info )
     IF( info /= 0 ) CALL errore( ' rdiaghg ', ' problems computing cholesky ', ABS( info ) )
#else
264
     CALL qe_pdpotrf( ss, nx, n, desc )
265
#endif
266 267 268
     !
  END IF
  !
269
  CALL stop_clock( 'rdiaghg:choldc' )
270 271 272
  !
  ! ... L is inverted ( s = L^-1 )
  !
273
  CALL start_clock( 'rdiaghg:inversion' )
274
  !
275
  IF( desc%active_node > 0 ) THEN
276
     !
277
#if defined(__SCALAPACK)
278 279
     ! 
     CALL sqr_dsetmat( 'U', n, zero, ss, size(ss,1), desc )
280 281 282 283 284

     CALL PDTRTRI( 'L', 'N', n, ss, 1, 1, desch, info )
     !
     IF( info /= 0 ) CALL errore( ' rdiaghg ', ' problems computing inverse ', ABS( info ) )
#else
285
     CALL qe_pdtrtri ( ss, nx, n, desc )
286
#endif
287 288 289
     !
  END IF
  !
290
  CALL stop_clock( 'rdiaghg:inversion' )
291 292 293
  !
  ! ... v = L^-1*H
  !
294
  CALL start_clock( 'rdiaghg:paragemm' )
295
  !
296
  IF( desc%active_node > 0 ) THEN
297 298 299 300 301 302 303
     !
     CALL sqr_mm_cannon( 'N', 'N', n, ONE, ss, nx, hh, nx, ZERO, v, nx, desc )
     !
  END IF
  !
  ! ... h = ( L^-1*H )*(L^-1)^T
  !
304
  IF( desc%active_node > 0 ) THEN
305 306 307 308 309
     !
     CALL sqr_mm_cannon( 'N', 'T', n, ONE, v, nx, ss, nx, ZERO, hh, nx, desc )
     !
  END IF
  !
310
  CALL stop_clock( 'rdiaghg:paragemm' )
311
  !
312
  IF ( desc%active_node > 0 ) THEN
313 314 315
     ! 
     !  Compute local dimension of the cyclically distributed matrix
     !
316
#if defined(__SCALAPACK)
ccavazzoni's avatar
ccavazzoni committed
317
     CALL pdsyevd_drv( .true., n, desc%nrcx, hh, SIZE(hh,1), e, ortho_cntx, ortho_comm )
318
#else
319
     CALL qe_pdsyevd( .true., n, desc, hh, SIZE(hh,1), e )
320
#endif
321 322 323 324 325
     !
  END IF
  !
  ! ... v = (L^T)^-1 v
  !
326
  CALL start_clock( 'rdiaghg:paragemm' )
327
  !
328
  IF ( desc%active_node > 0 ) THEN
329 330 331 332 333 334 335 336
     !
     CALL sqr_mm_cannon( 'T', 'N', n, ONE, ss, nx, hh, nx, ZERO, v, nx, desc )
     !
     DEALLOCATE( ss )
     DEALLOCATE( hh )
     !
  END IF
  !
337
  CALL mp_bcast( e, root, ortho_parent_comm )
338
  !
339
  CALL stop_clock( 'rdiaghg:paragemm' )
340
  !
341
  CALL stop_clock( 'rdiaghg' )
342
  !
343
  RETURN
344
  !
345
END SUBROUTINE prdiaghg