elpa.F90 12.7 KB
 Gregor Michalicek committed Jun 29, 2016 1 2 3 4 5 6 ``````!-------------------------------------------------------------------------------- ! 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. !-------------------------------------------------------------------------------- `````` Markus Betzinger committed Apr 26, 2016 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 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 ``````#ifdef CPP_INVERSION #define CPP_CHOLESKY cholesky_real #define CPP_invert_trm invert_trm_real #define CPP_solve_evp solve_evp_real #define CPP_solve_evp_2stage solve_evp_real_2stage #define CPP_REALCOMPLEX real #define CPP_transpose pdtran #define CPP_ONE 1.0 #define CPP_ZERO 0.0 #define CPP_mult mult_at_b_real #else #define CPP_CHOLESKY cholesky_complex #define CPP_invert_trm invert_trm_complex #define CPP_solve_evp solve_evp_complex #define CPP_solve_evp_2stage solve_evp_complex_2stage #define CPP_REALCOMPLEX complex #define CPP_transpose pztranc #define CPP_ONE cmplx(1.,0.) #define CPP_ZERO cmplx(0.,0.) #define CPP_mult mult_ah_b_complex #endif MODULE m_elpa PRIVATE !Module to call elpa library for parallel diagonalization !uses ssubredist1/2 for redistribution PUBLIC elpa CONTAINS SUBROUTINE elpa(m,n, SUB_COMM, a,b, z,eig,num) ! !---------------------------------------------------- !- Parallel eigensystem solver - driver routine based on chani; dw'12 ! ! m ........ actual (=leading) dimension of full a & b matrices ! must be problem size, as input a, b are one-dimensional ! and shall be redistributed to two-dimensional matrices ! actual (=leading) dimension of eigenvector z(,) ! n ........ number of columns of full (sub)matrix ( about n/np) ! SUB_COMM.. communicator for MPI ! a,b .... packed (sub)matrices, here expanded to non-packed ! z,eig .... eigenvectors and values, output ! num ...... number of ev's searched (and found) on this node ! On input, overall number of ev's searched, ! On output, local number of ev's found ! !---------------------------------------------------- ! USE m_juDFT USE elpa1 #ifdef CPP_ELPA2 USE elpa2 #endif IMPLICIT NONE INCLUDE 'mpif.h' INTEGER, INTENT (IN) :: m,n INTEGER, INTENT (IN) :: SUB_COMM INTEGER, INTENT (INOUT) :: num REAL, INTENT (OUT) :: eig(:) CPP_REALCOMPLEX, ALLOCATABLE, INTENT (INOUT) :: a(:),b(:) CPP_REALCOMPLEX, INTENT (OUT) :: z(:,:) ! !... Local variables ! INTEGER :: sc_desc(9) !SCALAPACK DESCRIPTOR INTEGER :: np,nb,myid,npcol,nprow,mycol,myrow,myrowssca INTEGER :: mycolssca,ictextblacs,mpi_comm_rows INTEGER :: mpi_comm_cols INTEGER :: num2 LOGICAL :: ok INTEGER :: err,ierr INTEGER :: i,k,n_col,n_row ! large distributed Matrices REAL,ALLOCATABLE :: eig2(:) CPP_REALCOMPLEX, ALLOCATABLE :: asca(:,:), bsca(:,:),tmp2(:,:) CPP_REALCOMPLEX, ALLOCATABLE :: eigvec(:,:) INTEGER, EXTERNAL :: numroc, indxl2g !SCALAPACK functions CALL priv_create_blacsgrid(sub_comm,m, np,nb,myid,npcol,nprow,mycol,myrow,myrowssca,mycolssca,& ictextblacs,sc_desc,mpi_comm_rows,mpi_comm_cols) ! Number of columns the local process gets in ScaLAPACK distribution ALLOCATE ( asca(myrowssca,mycolssca), stat=err ) IF (err/=0) CALL juDFT_error("allocate asca", calledby="chani_elpa") asca=0.0 ALLOCATE ( bsca(myrowssca,mycolssca), stat=err ) IF (err/=0) CALL juDFT_error("allocate bsca", calledby="chani_elpa") bsca=0.0 ! ! transfer to ScaLAPACK block cyclic 2d distribution with ! block-size nb ! !print *,"Before subredist" CALL subredist1(a,m,asca,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ) CALL subredist1(b,m,bsca,myrowssca,SUB_COMM, nprow, npcol, myid, ierr, nb ) ! for saving memory one might deallocate(a,b) !print *,"Before Allocate" ALLOCATE ( eig2(m), stat=err ) ! The eigenvalue array for ScaLAPACK IF (err.NE.0) CALL juDFT_error('Failed to allocated "eig2"', calledby ='elpa') ALLOCATE ( eigvec(myrowssca,mycolssca), stat=err ) ! Eigenvectors for ScaLAPACK IF (err.NE.0) CALL juDFT_error('Failed to allocated "eigvec"',calledby ='elpa') ALLOCATE ( tmp2(myrowssca,mycolssca), stat=err ) ! tmp_array IF (err.NE.0) CALL juDFT_error('Failed to allocated "tmp2"', calledby ='elpa') !print *,"Before elpa" !ELPA -start here ! Solve generalized problem ! ! 1. Calculate Cholesky factorization of Matrix B = U**T * U ! and invert triangular matrix U ! ! Please note: cholesky_complex/invert_trm_complex are not trimmed for speed. ! The only reason having them is that the Scalapack counterpart ! PDPOTRF very often fails on higher processor numbers for unknown reasons! #ifdef CPP_ELPA_NEW `````` Daniel Wortmann committed Jun 10, 2016 131 `````` CALL CPP_CHOLESKY (m,bsca,SIZE(bsca,1),nb, mpi_comm_rows,mpi_comm_cols,ok) `````` Markus Betzinger committed Apr 26, 2016 132 133 `````` CALL CPP_invert_trm(m, bsca, SIZE(bsca,1), nb, mpi_comm_rows, mpi_comm_cols,ok) #else `````` Daniel Wortmann committed Jun 10, 2016 134 `````` CALL CPP_CHOLESKY (m,bsca,SIZE(bsca,1),nb, mpi_comm_rows,mpi_comm_cols) `````` Markus Betzinger committed Apr 26, 2016 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 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 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 `````` CALL CPP_invert_trm(m, bsca, SIZE(bsca,1), nb, mpi_comm_rows, mpi_comm_cols) #endif ! 2. Calculate U**-T * A * U**-1 ! 2a. eigvec = U**-T * A ! A is only set in the upper half, solve_evp_real needs a full matrix ! Set lower half from upper half CALL CPP_transpose (m,m,CPP_ONE,asca,1,1,sc_desc, CPP_ZERO,eigvec,1,1,sc_desc) DO i=1,SIZE(asca,1) ! Get global column corresponding to i and number of local rows up to ! and including the diagonal, these are unchanged in A n_col = indxl2g(i, nb, mycol, 0, npcol) n_row = numroc (n_col, nb, myrow, 0, nprow) asca(n_row+1:myrowssca,i) = eigvec(n_row+1:myrowssca,i) ENDDO CALL CPP_mult ('U', 'L',m, m,bsca,myrowssca,asca,SIZE(asca,1),nb, mpi_comm_rows, mpi_comm_cols,eigvec,myrowssca) ! 2b. tmp2 = eigvec**T CALL CPP_transpose(m,m,CPP_ONE,eigvec,1,1,sc_desc,CPP_ZERO,tmp2,1,1,sc_desc) ! 2c. A = U**-T * tmp2 ( = U**-T * Aorig * U**-1 ) CALL CPP_mult ('U', 'U', m, m, bsca, SIZE(bsca,1), tmp2,& SIZE(tmp2,1),nb, mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1)) ! A is only set in the upper half, solve_evp_real needs a full matrix ! Set lower half from upper half CALL CPP_transpose(m,m,CPP_ONE,asca,1,1,sc_desc, CPP_ZERO,eigvec,1,1,sc_desc) DO i=1,SIZE(asca,1) ! Get global column corresponding to i and number of local rows up to ! and including the diagonal, these are unchanged in A n_col = indxl2g(i, nb, mycol, 0, npcol) n_row = numroc (n_col, nb, myrow, 0, nprow) asca(n_row+1:myrowssca,i) = eigvec(n_row+1:myrowssca,i) ENDDO ! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1 ! Eigenvectors go to eigvec num2=num #ifdef CPP_ELPA_NEW #ifdef CPP_ELPA2 err=CPP_solve_evp_2stage(m,num2,asca,SIZE(asca,1),& eig2,eigvec,SIZE(asca,1), nb, mpi_comm_rows, mpi_comm_cols,sub_comm) #else err=CPP_solve_evp(m, num2,asca,SIZE(asca,1),& eig2,eigvec, SIZE(asca,1), nb,mpi_comm_rows, mpi_comm_cols) #endif #else #ifdef CPP_ELPA2 CALL CPP_solve_evp_2stage(m,& num2,asca,SIZE(asca,1),eig2,eigvec,SIZE(asca,1), nb,mpi_comm_rows, mpi_comm_cols,sub_comm) #else CALL CPP_solve_evp(m, num2,& asca,SIZE(asca,1),eig2,eigvec, SIZE(asca,1), nb,mpi_comm_rows, mpi_comm_cols) #endif #endif ! 4. Backtransform eigenvectors: Z = U**-1 * eigvec ! mult_ah_b_complex needs the transpose of U**-1, thus tmp2 = (U**-1)**T CALL CPP_transpose(m,m,CPP_ONE,bsca,1,1,sc_desc,CPP_ZERO,tmp2,1,1,sc_desc) CALL CPP_mult ('L', 'N',m, num2, tmp2, SIZE(asca,1),& eigvec, SIZE(asca,1),nb,mpi_comm_rows, mpi_comm_cols, asca, SIZE(asca,1)) ! END of ELPA stuff CALL BLACS_GRIDEXIT(ictextblacs,ierr) CALL MPI_COMM_FREE(mpi_comm_rows,ierr) CALL MPI_COMM_FREE(mpi_comm_cols,ierr) !print *,"elpa done" ! ! Put those eigenvalues expected by chani to eig, i.e. for ! process i these are eigenvalues i+1, np+i+1, 2*np+i+1... ! Only num=num2/np eigenvalues per process ! num=FLOOR(REAL(num2)/np) IF (myid.LT.num2-(num2/np)*np) num=num+1 k=1 DO i=myid+1,num2,np eig(k)=eig2(i) k=k+1 ENDDO !write(*,"(a,i5,99f10.4)") "ELPA:",myid,eig ! ! Redistribute eigvec from ScaLAPACK distribution to each process ! having all eigenvectors corresponding to his eigenvalues as above ! CALL subredist2(z,m,num2,asca,myrowssca,SUB_COMM,nprow,npcol, myid,ierr,nb) ! DEALLOCATE ( asca ) DEALLOCATE ( bsca ) DEALLOCATE (eigvec, tmp2, eig2) END SUBROUTINE elpa SUBROUTINE priv_create_blacsgrid(mpi_subcom,m, np,nb,myid,npcol,nprow,mycol,& myrow,myrowssca,mycolssca, ictextblacs,sc_desc,mpi_comm_rows,mpi_comm_cols) USE m_juDFT USE elpa1 IMPLICIT NONE INCLUDE 'mpif.h' INTEGER,INTENT(IN) :: mpi_subcom INTEGER,INTENT(IN) :: m INTEGER,INTENT(OUT):: np,nb, npcol,nprow,myrowssca,mycolssca INTEGER,INTENT(OUT):: myrow,mycol,myid INTEGER,INTENT(OUT):: ictextblacs,sc_desc(:) INTEGER,INTENT(OUT):: mpi_comm_rows,mpi_comm_cols INTEGER :: iamblacs,npblacs INTEGER :: nprow2,npcol2,myrowblacs,mycolblacs INTEGER :: k,i,j INTEGER :: ierr INTEGER,ALLOCATABLE :: iblacsnums(:),ihelp(:),iusermap(:,:) EXTERNAL descinit, blacs_get EXTERNAL blacs_pinfo, blacs_gridinit !Determine rank and no of processors CALL MPI_COMM_RANK(mpi_subcom,myid,ierr) CALL MPI_COMM_SIZE(mpi_subcom,np,ierr) !print *,"priv_create_blacsgrid" ! determine block size ! nb = 20 IF (m.GT.2048) nb = 30 !2 IF (m.GT.8*2048) nb = 60 !4 ! compute processor grid, as square as possible ! If not square with more rows than columns distloop: DO j=INT(SQRT(REAL(np))),1,-1 IF ( (np/j) * j == np) THEN npcol = np/j nprow = j EXIT distloop ENDIF ENDDO distloop ALLOCATE(iblacsnums(np),ihelp(np),iusermap(nprow,npcol)) ! An nprow*npcol processor grid will be created ! Row and column index myrow, mycol of this processor in the grid ! and distribution of A and B in ScaLAPACK ! The local processor will get myrowssca rows and mycolssca columns ! of A and B ! myrow = myid/npcol ! my row number in the BLACS nprow*npcol grid mycol = myid -(myid/npcol)*npcol ! my column number in the BLACS nprow*npcol grid ! ! Now allocate Asca to put the elements of Achi or receivebuffer to ! myrowssca=(m-1)/(nb*nprow)*nb+ MIN(MAX(m-(m-1)/(nb*nprow)*nb*nprow-nb*myrow,0),nb) ! Number of rows the local process gets in ScaLAPACK distribution mycolssca=(m-1)/(nb*npcol)*nb+ MIN(MAX(m-(m-1)/(nb*npcol)*nb*npcol-nb*mycol,0),nb) !Get BLACS ranks for all MPI ranks CALL BLACS_PINFO(iamblacs,npblacs) ! iamblacs = local process rank (e.g. myid) ! npblacs = number of available processes iblacsnums=-2 ihelp=-2 ihelp(myid+1)=iamblacs ! Get the Blacs id corresponding to the MPI id !print *,"ALLREDUCE:",mpi_subcom CALL MPI_ALLREDUCE(ihelp, iblacsnums, np,MPI_INTEGER,MPI_MAX,mpi_subcom,ierr) IF (ierr.NE.0)CALL juDFT_error('Error in allreduce for BLACS nums' ,calledby='elpa') ! iblacsnums(i) is the BLACS-process number of MPI-process i-1 k = 1 DO i = 1, nprow DO j = 1, npcol iusermap(i,j) = iblacsnums(k) k = k + 1 ENDDO ENDDO !Get the Blacs default context CALL BLACS_GET(0,0,ictextblacs) ! Create the Grid CALL BLACS_GRIDMAP(ictextblacs,iusermap,nprow,nprow,npcol) ! Now control, whether the BLACS grid is the one we wanted CALL BLACS_GRIDINFO(ictextblacs, nprow2,npcol2,myrowblacs,mycolblacs) IF (nprow2 /= nprow) THEN WRITE(6,*) 'Wrong number of rows in BLACS grid' WRITE(6,*) 'nprow=',nprow,' nprow2=',nprow2 CALL juDFT_error('Wrong number of rows in BLACS grid',calledby= 'elpa') ENDIF IF (npcol2 /= npcol) THEN WRITE(6,*) 'Wrong number of columns in BLACS grid' WRITE(6,*) 'npcol=',npcol,' npcol2=',npcol2 CALL juDFT_error('Wrong number of columns in BLACS grid', calledby ='elpa') ENDIF !Create communicators for ELPA !print *,"creating ELPA comms" CALL get_elpa_row_col_comms(mpi_subcom, myrowblacs, mycolblacs,mpi_comm_rows, mpi_comm_cols) !print *,"creating ELPA comms -- done" !Create the descriptors CALL descinit(sc_desc,m,m,nb,nb,0,0,ictextblacs,myrowssca,ierr) IF (ierr /=0 ) CALL juDFT_error('descinit1 failed',calledby='elpa') END SUBROUTINE priv_create_blacsgrid END MODULE m_elpa``````