scalapack.F90 10.3 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_scalapack
CONTAINS
  SUBROUTINE scalapack(hmat,smat,ne,eig,ev)
    !
    !----------------------------------------------------
    !- Parallel eigensystem solver - driver routine; gb99
    !  Uses the SCALAPACK for the actual diagonalization
    !
Uliana Alekseeva's avatar
Uliana Alekseeva committed
15 16 17
    ! hmat ..... Hamiltonian matrix
    ! smat ..... overlap matrix
    ! ne ....... number of ev's searched (and found) on this node
Daniel Wortmann's avatar
Daniel Wortmann committed
18 19
    !            On input, overall number of ev's searched,
    !            On output, local number of ev's found
20 21
    ! eig ...... all eigenvalues, output
    ! ev ....... local eigenvectors, output
Daniel Wortmann's avatar
Daniel Wortmann committed
22 23 24 25 26
    !
    !----------------------------------------------------
    !
#include"cpp_double.h"
    USE m_juDFT
27
    USE m_constants
Daniel Wortmann's avatar
Daniel Wortmann committed
28
    USE m_types_mpimat
29
    USE m_types_mat
Daniel Wortmann's avatar
Daniel Wortmann committed
30 31 32 33 34 35 36 37 38 39 40 41 42
    IMPLICIT NONE
    CLASS(t_mat),INTENT(INOUT)    :: hmat,smat
    CLASS(t_mat),ALLOCATABLE,INTENT(OUT)::ev
    REAL,INTENT(out)              :: eig(:)
    INTEGER,INTENT(INOUT)         :: ne
    
    
#ifdef CPP_SCALAPACK
#ifdef CPP_MPI    
    INCLUDE 'mpif.h'
#endif
    !...  Local variables
    !
Uliana Alekseeva's avatar
Uliana Alekseeva committed
43 44
    INTEGER i , ierr, err
    INTEGER, ALLOCATABLE :: iwork(:)
Daniel Wortmann's avatar
Daniel Wortmann committed
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
    REAL,    ALLOCATABLE :: rwork(:)
    INTEGER              :: lrwork

    !
    !  ScaLAPACK things
    CHARACTER (len=1)    :: uplo
    INTEGER              :: num,num1,num2,liwork,lwork2,np0,mq0,np,myid
    INTEGER              :: iceil, numroc, nn, nb
    INTEGER, ALLOCATABLE :: ifail(:), iclustr(:)
    REAL                 :: abstol,orfac=1.E-4,CPP_LAPACK_slamch
    REAL,ALLOCATABLE     :: eig2(:), gap(:)
    REAL,    ALLOCATABLE :: work2_r(:)
    COMPLEX, ALLOCATABLE :: work2_c(:)

    TYPE(t_mpimat):: ev_dist
    
    EXTERNAL iceil, numroc
    EXTERNAL CPP_LAPACK_slamch
    
    
    SELECT TYPE(hmat)
    TYPE IS (t_mpimat)
    SELECT TYPE(smat)
    TYPE IS (t_mpimat)

    ALLOCATE(eig2(hmat%global_size1))


73 74
    CALL MPI_COMM_RANK(hmat%blacsdata%mpi_com,myid,ierr)
    CALL MPI_COMM_SIZE(hmat%blacsdata%mpi_com,np,ierr)
Daniel Wortmann's avatar
Daniel Wortmann committed
75 76 77 78 79

    num=ne !no of states solved for
    
    abstol=2.0*CPP_LAPACK_slamch('S') ! PDLAMCH gave an error on ZAMpano

80
    CALL ev_dist%init(hmat)
Daniel Wortmann's avatar
Daniel Wortmann committed
81

82 83
    !smat%blacs_desc(2)    = hmat%blacs_desc(2)
    !ev_dist%blacs_desc(2) = hmat%blacs_desc(2)
84 85
    !smat%blacs_desc=hmat%blacs_desc
    !ev_dist%blacs_desc=hmat%blacs_desc
86 87
    
    
88 89
    nb=hmat%blacsdata%blacs_desc(5)! Blocking factor
    IF (nb.NE.hmat%blacsdata%blacs_desc(6)) CALL judft_error("Different block sizes for rows/columns not supported")
Daniel Wortmann's avatar
Daniel Wortmann committed
90 91 92
  
    !
    nn=MAX(MAX(hmat%global_size1,nb),2)
93 94
    np0=numroc(nn,nb,0,0,hmat%blacsdata%nprow)
    mq0=numroc(MAX(MAX(ne,nb),2),nb,0,0,hmat%blacsdata%npcol)
Daniel Wortmann's avatar
Daniel Wortmann committed
95
    IF (hmat%l_real) THEN
96
       lwork2=5*hmat%global_size1+MAX(5*nn,np0*mq0+2*nb*nb)+ iceil(ne,hmat%blacsdata%nprow*hmat%blacsdata%npcol)*nn
Daniel Wortmann's avatar
Daniel Wortmann committed
97 98 99 100 101 102 103 104 105 106
       ALLOCATE ( work2_r(lwork2+10*hmat%global_size1), stat=err ) ! Allocate more in case of clusters
    ELSE
       lwork2=hmat%global_size1+MAX(nb*(np0+1),3)
       ALLOCATE ( work2_c(lwork2), stat=err )
    ENDIF
    IF (err.NE.0) THEN 
       WRITE (*,*) 'work2  :',err,lwork2
       CALL juDFT_error('Failed to allocated "work2"', calledby ='chani')
    ENDIF
    
107
    liwork=6*MAX(MAX(hmat%global_size1,hmat%blacsdata%nprow*hmat%blacsdata%npcol+1),4)
Daniel Wortmann's avatar
Daniel Wortmann committed
108 109 110 111 112 113 114 115 116 117
    ALLOCATE ( iwork(liwork), stat=err )
    IF (err.NE.0) THEN
       WRITE (*,*) 'iwork  :',err,liwork
       CALL juDFT_error('Failed to allocated "iwork"', calledby ='chani')
    ENDIF
    ALLOCATE ( ifail(hmat%global_size1), stat=err )
    IF (err.NE.0) THEN
       WRITE (*,*) 'ifail  :',err,hmat%global_size1
       CALL juDFT_error('Failed to allocated "ifail"', calledby ='chani')
    ENDIF
118
    ALLOCATE ( iclustr(2*hmat%blacsdata%nprow*hmat%blacsdata%npcol), stat=err )
Daniel Wortmann's avatar
Daniel Wortmann committed
119
    IF (err.NE.0) THEN
120
       WRITE (*,*) 'iclustr:',err,2*hmat%blacsdata%nprow*hmat%blacsdata%npcol
Daniel Wortmann's avatar
Daniel Wortmann committed
121 122
       CALL juDFT_error('Failed to allocated "iclustr"', calledby ='chani')
    ENDIF
123
    ALLOCATE ( gap(hmat%blacsdata%nprow*hmat%blacsdata%npcol), stat=err )
Daniel Wortmann's avatar
Daniel Wortmann committed
124
    IF (err.NE.0) THEN
125
       WRITE (*,*) 'gap    :',err,hmat%blacsdata%nprow*hmat%blacsdata%npcol
Daniel Wortmann's avatar
Daniel Wortmann committed
126 127 128 129 130 131 132 133
       CALL juDFT_error('Failed to allocated "gap"', calledby ='chani')
    ENDIF
    !
    !     Compute size of workspace
    !
    IF (hmat%l_real) THEN
       uplo='U'
       CALL CPP_LAPACK_pdsygvx(1,'V','I','U',hmat%global_size1,hmat%data_r,1,1,&
134
            hmat%blacsdata%blacs_desc,smat%data_r,1,1,smat%blacsdata%blacs_desc,&
Daniel Wortmann's avatar
Daniel Wortmann committed
135
            0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,ev_dist%data_r,1,1,&
136
            ev_dist%blacsdata%blacs_desc,work2_r,-1,iwork,-1,ifail,iclustr, gap,ierr)
Daniel Wortmann's avatar
Daniel Wortmann committed
137 138 139 140 141 142 143 144 145 146
       IF ( work2_r(1).GT.lwork2) THEN
          lwork2 = work2_r(1)
          DEALLOCATE (work2_r)
          ALLOCATE ( work2_r(lwork2+20*hmat%global_size1), stat=err ) ! Allocate even more in case of clusters
          IF (err.NE.0) THEN
             WRITE (*,*) 'work2  :',err,lwork2
             CALL juDFT_error('Failed to allocated "work2"', calledby ='chani')
          ENDIF
       ENDIF
    ELSE
147
       lrwork=4*hmat%global_size1+MAX(5*nn,np0*mq0)+ iceil(ne,hmat%blacsdata%nprow*hmat%blacsdata%npcol)*nn
Daniel Wortmann's avatar
Daniel Wortmann committed
148 149 150 151 152 153 154 155
       ! Allocate more in case of clusters
       ALLOCATE(rwork(lrwork+10*hmat%global_size1), stat=ierr)
       IF (err /= 0) THEN
          WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed'
          CALL juDFT_error('Failed to allocated "rwork"', calledby ='chani')
       ENDIF
       
       CALL CPP_LAPACK_pzhegvx(1,'V','I','U',hmat%global_size1,hmat%data_c,1,1,&
156
            hmat%blacsdata%blacs_desc,smat%data_c,1,1, smat%blacsdata%blacs_desc,&
Daniel Wortmann's avatar
Daniel Wortmann committed
157
            0.0,1.0,1,num,abstol,num1,num2,eig2,orfac,ev_dist%data_c,1,1,&
158
            ev_dist%blacsdata%blacs_desc,work2_c,-1,rwork,-1,iwork,-1,ifail,iclustr,&
Daniel Wortmann's avatar
Daniel Wortmann committed
159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
            gap,ierr)
       IF (ABS(work2_c(1)).GT.lwork2) THEN
          lwork2=work2_c(1)
          DEALLOCATE (work2_c)
          ALLOCATE (work2_c(lwork2), stat=err)
          IF (err /= 0) THEN
             WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed:',lwork2
             CALL juDFT_error('Failed to allocated "work2"', calledby ='chani')
          ENDIF
       ENDIF
       IF (rwork(1).GT.lrwork) THEN
          lrwork=rwork(1)
          DEALLOCATE(rwork)
          ! Allocate even more in case of clusters
          ALLOCATE (rwork(lrwork+20*hmat%global_size1), stat=err)
          IF (err /= 0) THEN
             WRITE (*,*) 'ERROR: chani.F: Allocating rwork failed: ', lrwork+20*hmat%global_size1
             CALL juDFT_error('Failed to allocated "rwork"', calledby ='chani')
          ENDIF
       ENDIF
    endif
    IF (iwork(1) .GT. liwork) THEN
       liwork = iwork(1)
       DEALLOCATE (iwork)
       ALLOCATE (iwork(liwork), stat=err)
       IF (err /= 0) THEN
          WRITE (*,*) 'ERROR: chani.F: Allocating iwork failed: ',liwork
          CALL juDFT_error('Failed to allocated "iwork"', calledby ='chani')
       ENDIF
    ENDIF
    !
    !     Now solve generalized eigenvalue problem
    !
    CALL timestart("SCALAPACK call")
    if (hmat%l_real) THEN
194
       CALL CPP_LAPACK_pdsygvx(1,'V','I','U',hmat%global_size1,hmat%data_r,1,1,hmat%blacsdata%blacs_desc,smat%data_r,1,1, smat%blacsdata%blacs_desc,&
Daniel Wortmann's avatar
Daniel Wortmann committed
195
            1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,ev_dist%data_r,1,1,&
196
            ev_dist%blacsdata%blacs_desc,work2_r,lwork2,iwork,liwork,ifail,iclustr,&
Daniel Wortmann's avatar
Daniel Wortmann committed
197 198
            gap,ierr)
    else
199
       CALL CPP_LAPACK_pzhegvx(1,'V','I','U',hmat%global_size1,hmat%data_c,1,1,hmat%blacsdata%blacs_desc,smat%data_c,1,1, smat%blacsdata%blacs_desc,&
Daniel Wortmann's avatar
Daniel Wortmann committed
200
            1.0,1.0,1,num,abstol,num1,num2,eig2,orfac,ev_dist%data_c,1,1,&
201
            ev_dist%blacsdata%blacs_desc,work2_c,lwork2,rwork,lrwork,iwork,liwork,&
Daniel Wortmann's avatar
Daniel Wortmann committed
202 203 204 205 206
            ifail,iclustr,gap,ierr)
       DEALLOCATE(rwork)
    endif
    CALL timestop("SCALAPACK call")
    IF (ierr .NE. 0) THEN
207 208
       !IF (ierr /= 2) WRITE (oUnit,*) myid,' error in pzhegvx/pdsygvx, ierr=',ierr
       !IF (ierr <= 0) WRITE (oUnit,*) myid,' illegal input argument'
Daniel Wortmann's avatar
Daniel Wortmann committed
209
       IF (MOD(ierr,2) /= 0) THEN
210
          !WRITE (oUnit,*) myid,'some eigenvectors failed to converge'
Daniel Wortmann's avatar
Daniel Wortmann committed
211 212
          eigs: DO i = 1, ne
             IF (ifail(i) /= 0) THEN
213
                !WRITE (oUnit,*) myid,' eigenvector',ifail(i), 'failed to converge'
Daniel Wortmann's avatar
Daniel Wortmann committed
214 215 216 217
             ELSE
                EXIT eigs
             ENDIF
          ENDDO eigs
218
          !CALL CPP_flush(oUnit)
Daniel Wortmann's avatar
Daniel Wortmann committed
219 220
       ENDIF
       IF (MOD(ierr/4,2).NE.0) THEN
221 222
          !WRITE(oUnit,*) myid,' only',num2,' eigenvectors converged'
          !CALL CPP_flush(oUnit)
Daniel Wortmann's avatar
Daniel Wortmann committed
223 224
       ENDIF
       IF (MOD(ierr/8,2).NE.0) THEN
225
          !WRITE(oUnit,*) myid,' PDSTEBZ failed to compute eigenvalues'
Daniel Wortmann's avatar
Daniel Wortmann committed
226 227 228
          CALL judft_warn("SCALAPACK failed to solve eigenvalue problem",calledby="scalapack.f90")
       ENDIF
       IF (MOD(ierr/16,2).NE.0) THEN
229
          !WRITE(oUnit,*) myid,' B was not positive definite, Cholesky failed at',ifail(1)
Daniel Wortmann's avatar
Daniel Wortmann committed
230 231 232 233 234
          CALL judft_warn("SCALAPACK failed: B was not positive definite",calledby="scalapack.f90")
       ENDIF
    ENDIF
    IF (num2 < num1) THEN
       !IF (myid ==0) THEN
235 236 237 238
          WRITE(oUnit,*) 'Not all eigenvalues wanted are found'
          WRITE(oUnit,*) 'number of eigenvalues/vectors wanted',num1
          WRITE(oUnit,*) 'number of eigenvalues/vectors found',num2
          !CALL CPP_flush(oUnit)
Daniel Wortmann's avatar
Daniel Wortmann committed
239 240 241
       !ENDIF
    ENDIF
    !
242 243 244 245 246 247 248 249
    !     Each process has all eigenvalues in output
    eig(:num2) = eig2(:num2)    
    DEALLOCATE(eig2)
    !
    !
    !     Redistribute eigenvectors  from ScaLAPACK distribution to each process, i.e. for
    !     process i these are eigenvectors i+1, np+i+1, 2*np+i+1...
    !     Only num=num2/np eigenvectors per process
Daniel Wortmann's avatar
Daniel Wortmann committed
250 251 252 253 254 255
    !
    num=FLOOR(REAL(num2)/np)
    IF (myid.LT.num2-(num2/np)*np) num=num+1
    ne=0
    DO i=myid+1,num2,np
       ne=ne+1
256
       !eig(ne)=eig2(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
257 258
    ENDDO
    ALLOCATE(t_mpimat::ev)
259
    CALL ev%init(ev_dist%l_real,ev_dist%global_size1,ev_dist%global_size1,ev_dist%blacsdata%mpi_com,.FALSE.)
Daniel Wortmann's avatar
Daniel Wortmann committed
260 261 262 263 264 265 266 267 268 269
    CALL ev%copy(ev_dist,1,1)
    CLASS DEFAULT
      call judft_error("Wrong type (1) in scalapack")
    END SELECT
    CLASS DEFAULT
      call judft_error("Wrong type (2) in scalapack")
    END SELECT
#endif
  END SUBROUTINE scalapack
END MODULE m_scalapack