2
! Copyright (C) 1996-2016 The SIESTA group
3
! This file is distributed under the terms of the
4
! GNU General Public License: see COPYING in the top directory
5
! or http://www.gnu.org/copyleft/gpl.txt .
6
! See Docs/Contributors.txt for a list of contributors.
11
#ifndef SIESTA__NO_MRRR
15
subroutine rdiag_mrrr(H,S,n,nm,nml,w,Z,neigvec,ierror)
16
! ***************************************************************************
17
! Subroutine to find a subset of eigenvalues and eigenvectors of the
18
! real general eigenvalue problem H z = w S z, with H and S
19
! real symmetric matrices. It uses the MRRR algorithm.
21
! Written by A. Garcia (Jan 2014), based on rdiag.
23
! ************************** INPUT ******************************************
24
! real*8 H(nml,nm) : Symmetric H matrix
25
! real*8 S(nml,nm) : Symmetric S matrix
26
! integer n : Order of the generalized system
27
! integer nm : Right hand dimension of H and S matrices
28
! integer nml : Left hand dimension of H and S matrices
29
! which is greater than or equal to nm
30
! integer neigvec : No. of eigenvectors to calculate
31
! ************************** OUTPUT *****************************************
32
! real*8 w(nml) : Eigenvalues
33
! real*8 Z(nml,nm) : Eigenvectors
34
! integer ierror : Flag indicating success code for routine
40
use parallel, only : Node, Nodes, BlockSize
42
use mpi_siesta, only : mpi_comm_world
50
integer, intent(out) :: ierror
51
integer, intent(in) :: n
52
integer, intent(in) :: neigvec
53
integer, intent(in) :: nm
54
integer, intent(in) :: nml
55
real(dp), intent(in) :: H(nml,nm)
56
real(dp), intent(in) :: S(nml,nm)
58
! It should be possible to dimension Z to use the
59
! actual number of eigenvectors needed...
60
real(dp), intent(out) :: w(nml)
61
real(dp), intent(out) :: Z(nml,nm)
65
integer :: MPIerror, mpi_comm_rows, mpi_comm_cols
73
! Additional variables for a 2D proc grid
74
integer :: np_rows, np_cols
75
integer :: my_prow, my_pcol
77
integer, dimension(9) :: desc_h2d
79
integer :: neigok, nz, vl, vu
80
character(len=256) :: msg
83
real(dp), pointer, dimension(:,:) :: h2d =>null()
84
real(dp), pointer, dimension(:,:) :: s2d =>null()
85
real(dp), pointer, dimension(:,:) :: z2d =>null()
86
integer :: nh_rows, nh_cols
88
integer :: info, i, n_col, n_row
89
integer :: lwork, liwork
90
real(dp), pointer :: work(:) => null()
91
integer , pointer :: iwork(:) => null()
93
real(dp), parameter :: zero = 0.0_dp
94
real(dp), parameter :: one = 1.0_dp
96
integer, external :: numroc, indxl2g
98
!***********************************************************************
100
!***********************************************************************
102
! Initialise error flag
107
w(1) = H(1,1) / S(1,1)
109
Z(1,1) = one / sqrt(S(1,1))
114
call timer('rdiag_mrrr',1)
117
! Get Blacs context and initialise Blacs grid for main grid
119
! First, declare explicitly a context for the original
120
! distribution, which is block-cyclic over the columns,
121
! since our matrices are dimensioned as H(n,no_l) in diagg,
122
! i.e., each processor contains all the rows but just no_l
126
call blacs_get( -1, 0, ictxt )
127
! The design of the BLACS is weird.
128
! The above call is equivalent to setting:
129
ictxt = mpi_comm_world ! or whatever comm we are using
130
! The following call will OVERWRITE ictxt with a new context
131
! "C" means "column-major" ordering, which is not surprising,
132
! since we have np_rows_1d = 1
133
call blacs_gridinit( ictxt, 'C', np_rows_1d, np_cols_1d )
136
! Set up blacs descriptors for parallel case
138
! This is wrong in the BSs
139
! since the effective row blocksize is n. It should be
140
! n, n, n , BlockSize
141
call descinit( desch, n, n, BlockSize, BlockSize, 0, 0,
143
if (info.ne.0) BlacsOK = .false.
144
if (.not.BlacsOK) then
145
call die('ERROR : Blacs setup has failed in rdiag!')
148
! Setup secondary grid for 2D distribution
150
! Selection of number of processor rows/columns
151
! We try to set up the grid square-like, i.e. start the search for possible
152
! divisors of Nodes with a number next to the square root of Nodes
153
! and decrement it until a divisor is found.
155
do np_rows = NINT(SQRT(REAL(Nodes))),2,-1
156
if(mod(Nodes,np_rows) == 0 ) exit
158
! at the end of the above loop, Nodes is always divisible by np_rows
160
np_cols = Nodes/np_rows
162
!if (Node == 0) print *, "np_rows, col: ", np_rows, np_cols
164
!call blacs_get(ictxt, 10, i2d_ctxt)
165
i2d_ctxt = mpi_comm_world ! or whatever comm we are using
166
! after the above call, i2d_ctxt will contain the communicator
167
! handle we used to initialize ictxt
168
! Now we set up the process grid in row-major, or "natural" order.
171
call blacs_gridinit(i2d_ctxt, 'R', np_rows, np_cols)
173
! We need this call to know who we are (my_prow, my_pcol)
174
call blacs_gridinfo(i2d_ctxt, np_rows, np_cols,
176
!if (Node == 0) print *, "ictxt, i2_ctxt:", ictxt, i2d_ctxt
178
! Enquire size of local part of 2D matrices
179
! since we are going to change and in principle we do not know.
181
! We might want to change the value of BlockSize...
183
nh_rows = numroc(n, BlockSize, my_prow, 0, np_rows)
184
nh_cols = numroc(n, BlockSize, my_pcol, 0, np_cols)
186
! Set up blacs descriptors for 2D case
187
! It is only necessary to give the leading (row) dimension of the array.
188
call descinit(desc_h2d, n, n, Blocksize, BlockSize, 0, 0,
189
. i2d_ctxt, nh_rows, info)
190
if (info.ne.0) BlacsOK = .false.
191
if (.not.BlacsOK) then
192
call die('ERROR : Blacs setup has failed in rdiag!')
195
! Set up workspace arrays for 2D versions of 1D arrays
196
call re_alloc( h2d, 1, nh_rows, 1, nh_cols, name='h2d')
197
call re_alloc( s2d, 1, nh_rows, 1, nh_cols, name='s2d')
198
call re_alloc( z2d, 1, nh_rows, 1, nh_cols, name='z2d')
200
! Copy arrays to new 2D distribution
201
call pdgemr2d(n, n, H, 1, 1, desch, h2d, 1, 1, desc_h2d, ictxt)
202
call pdgemr2d(n, n, S, 1, 1, desch, s2d, 1, 1, desc_h2d, ictxt)
204
np0 = numroc(n,BlockSize,0,0,np_rows)
205
! lwork is used also in the forward transform.
206
! It will be reset below
207
lwork = max(BlockSize*(np0+1),3*BlockSize)
208
liwork = 3 ! for now: see query below
210
! Allocate workspace arrays
211
call re_alloc( work, 1,lwork, name='work' )
212
call re_alloc( iwork, 1,liwork, name='iwork' )
214
!*************************************************************************
215
! Transform problem to standard eigenvalue problem
216
! by first factorizing the overlap matrix by Cholesky
217
!*************************************************************************
218
call timer('rdiag1',1)
220
call pdpotrf('U',n,s2d,1,1,desc_h2d,info)
222
call die('Error in Cholesky factorisation in rdiag')
224
call timer('rdiag1',2)
229
call timer('rdiag2',1)
231
call pdsyngst(1,'U',n,h2d,1,1,desc_h2d,s2d,1,1,
232
. desc_h2d,dscale,work,lwork,info)
234
call die('Error in forward transformation in rdiag')
236
call timer('rdiag2',2)
238
! Solve standard eigenvalue problem
239
call timer('rdiag3',1)
241
! Find the optimal work-space sizes
244
call pdsyevr('V','I','U',n,h2d,1,1,desc_h2d,vl,vu,1,
245
. neigvec,neigok,nz,w,z2d,1,1,
246
. desc_h2d,work,lwork,iwork,liwork,info)
248
write(msg,"(a,i0)") "Error in pdsyevr work estim. Info: ",
252
lwork = nint(work(1))
254
call re_alloc(work, 1, lwork)
255
call re_alloc(iwork, 1, liwork)
257
call pdsyevr('V','I','U',n,h2d,1,1,desc_h2d,vl,vu,1,
258
. neigvec,neigok,nz,w,z2d,1,1,
259
. desc_h2d,work,lwork,iwork,liwork,info)
262
write(msg,"(a,i0)") "Error in pdsyevr. Info: ", info
266
if ((neigok /= neigvec)) then
267
write(6,"(a,2i7)") "*** Warning: eigenpair mismatch:"
269
$ "Eigenvalues requested and found (w,z): ",
270
$ neigvec , neigok, nz
271
write(0,"(a,2i7)") "*** Warning: eigenpair mismatch:"
273
$ "Eigenvalues requested and found (w,z): ",
274
$ neigvec , neigok, nz
278
call timer('rdiag3',2)
280
!********************************************************************
281
! Back transformation *
282
!********************************************************************
283
call timer('rdiag4',1)
285
call pdtrsm('Left','U','N','Non-unit',n,neigvec,one,
286
. s2d,1,1,desc_h2d,z2d,1,1,desc_h2d)
287
! Return to 1-D block-cyclic distribution
288
call pdgemr2d(n,n,z2d,1,1,desc_h2d,Z,1,1,desch,ictxt)
289
if (dscale.ne.one) then
290
call dscal(n,dscale,w,1)
294
call die('Error in back transformation in rdiag')
296
call timer('rdiag4',2)
298
CALL BLACS_GRIDEXIT( ICTXT )
299
CALL BLACS_GRIDEXIT( I2D_CTXT )
301
! Deallocate workspace arrays
302
call de_alloc( work, name='work' )
303
call de_alloc( iwork, name='iwork' )
306
call de_alloc( h2d, name='h2d' )
307
call de_alloc( s2d, name='s2d' )
308
call de_alloc( z2d, name='z2d' )
310
call timer('rdiag_mrrr',2)
313
end subroutine rdiag_mrrr
315
! Some compilers do not allow empty modules
317
subroutine rdiag_mrrr_dummy()
320
end module m_rdiag_mrrr
added
removed
Src/rdiag_mrrr.F
