|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00015 MODULE cp_fm_diag 00016 00017 USE cp_blacs_calls, ONLY: cp_blacs_gridexit,& 00018 cp_blacs_gridinit 00019 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 00020 cp_blacs_env_release 00021 USE cp_fm_basic_linalg, ONLY: cp_fm_syrk,& 00022 cp_fm_upper_to_full 00023 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00024 cp_fm_struct_release,& 00025 cp_fm_struct_type 00026 USE cp_fm_types, ONLY: cp_fm_create,& 00027 cp_fm_release,& 00028 cp_fm_to_fm,& 00029 cp_fm_type 00030 USE cp_para_env, ONLY: cp_para_env_create,& 00031 cp_para_env_release 00032 USE cp_para_types, ONLY: cp_blacs_env_type,& 00033 cp_para_env_type 00034 #if defined (__ELPA) 00035 USE ELPA2 00036 USE ELPA1 00037 #endif 00038 USE kinds, ONLY: dp,& 00039 dp_size,& 00040 int_size 00041 USE message_passing, ONLY: mp_bcast,& 00042 mp_comm_free,& 00043 mp_comm_split,& 00044 mp_comm_split_direct,& 00045 mp_sync 00046 USE termination, ONLY: stop_program 00047 USE timings, ONLY: timeset,& 00048 timestop 00049 #include "cp_common_uses.h" 00050 00051 IMPLICIT NONE 00052 00053 PRIVATE 00054 00055 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_diag' 00056 00057 ! Public subroutines 00058 00059 PUBLIC :: cp_fm_block_jacobi,& 00060 cp_fm_power, & 00061 cp_fm_syevd, & 00062 cp_fm_syevr, & 00063 cp_fm_syevx, & 00064 cp_fm_elpa, choose_eigv_solver 00065 00066 CONTAINS 00067 00068 00069 ! ***************************************************************************** 00076 SUBROUTINE choose_eigv_solver(matrix,eigenvectors,eigenvalues, info, error) 00077 00078 TYPE(cp_fm_type), POINTER :: matrix, eigenvectors 00079 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues 00080 INTEGER, INTENT(OUT), OPTIONAL :: info 00081 TYPE(cp_error_type), INTENT(inout) :: error 00082 00083 CHARACTER(len=*), PARAMETER :: routineN = 'choose_eigv_solver', 00084 routineP = moduleN//':'//routineN 00085 00086 INTEGER :: nmo 00087 LOGICAL :: use_elpa, use_SL2 00088 00089 IF(PRESENT(info)) info = 0 00090 00091 #if defined (__ELPA) 00092 use_elpa =.TRUE. 00093 #else 00094 use_elpa = .FALSE. 00095 #endif 00096 00097 #if defined (__SCALAPACK2) 00098 use_SL2 = .TRUE. 00099 #else 00100 use_SL2 = .FALSE. 00101 #endif 00102 00103 nmo=SIZE(eigenvalues,1) 00104 00105 IF(use_elpa) THEN 00106 CALL cp_fm_elpa(matrix,eigenvectors,eigenvalues,error=error) 00107 ELSEIF(use_SL2) THEN 00108 CALL cp_fm_syevr(matrix,eigenvectors,eigenvalues,1,nmo,error=error) 00109 ELSE 00110 CALL cp_fm_syevd(matrix,eigenvectors,eigenvalues,info=info,error=error) 00111 END IF 00112 00113 00114 END SUBROUTINE choose_eigv_solver 00115 00116 ! ***************************************************************************** 00125 SUBROUTINE cp_fm_syevd(matrix,eigenvectors,eigenvalues,info,error) 00126 00127 TYPE(cp_fm_type), POINTER :: matrix, eigenvectors 00128 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues 00129 INTEGER, INTENT(OUT), OPTIONAL :: info 00130 TYPE(cp_error_type), INTENT(inout) :: error 00131 00132 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_syevd', 00133 routineP = moduleN//':'//routineN 00134 00135 INTEGER :: fake_descriptor(9), handle, istat, liwork, lwork, mepos_old, 00136 myinfo, n, ngroups, nmo, num_pe_new, num_pe_old, subgroup 00137 INTEGER, DIMENSION(:), POINTER :: group_distribution, 00138 group_partition, iwork 00139 LOGICAL :: failure 00140 REAL(KIND=dp) :: fake_local_data(1,1) 00141 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig 00142 REAL(KIND=dp), DIMENSION(:), POINTER :: work 00143 REAL(KIND=dp), DIMENSION(:, :), POINTER :: m, v 00144 TYPE(cp_blacs_env_type), POINTER :: blacs_env_new 00145 TYPE(cp_fm_struct_type), POINTER :: fm_struct_new 00146 TYPE(cp_fm_type), POINTER :: eigenvectors_new, matrix_new 00147 TYPE(cp_para_env_type), POINTER :: para_env, para_env_new 00148 00149 ! ------------------------------------------------------------------------- 00150 00151 CALL timeset(routineN,handle) 00152 00153 failure = .FALSE. 00154 IF(PRESENT(info)) info = 0 00155 00156 n = matrix%matrix_struct%nrow_global 00157 ALLOCATE(eig(n), STAT=istat) 00158 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00159 00160 00161 #if defined(__SCALAPACK) 00162 00163 ! first figure out the optimal number of cpus 00164 ! this is pure heuristics, based on rosa timings 00165 ! that demonstrate that timings go up sharply if too many tasks are used 00166 ! we take a multiple of 4, and approximately n/60 00167 num_pe_old=matrix%matrix_struct%para_env%num_pe 00168 num_pe_new=((n+60*4-1)/(60*4))*4 00169 para_env=>matrix%matrix_struct%para_env 00170 mepos_old=para_env%mepos 00171 00172 ! if the optimal is smaller than num_pe, we will redistribute the input matrix 00173 ! it seems sensible to refactor redistributing matrices as this might be useful elsewhere 00174 IF (num_pe_new<num_pe_old) THEN 00175 00176 00177 ! split comm, the first num_pe_new tasks will do the work 00178 ALLOCATE(group_distribution(0:num_pe_old-1)) 00179 ALLOCATE(group_partition(0:1)) 00180 group_partition=(/num_pe_new,num_pe_old-num_pe_new/) 00181 CALL mp_comm_split(comm=para_env%group,sub_comm=subgroup,& 00182 ngroups=ngroups,group_distribution=group_distribution,& 00183 n_subgroups=2, group_partition=group_partition) 00184 00185 IF (group_distribution(mepos_old)==0) THEN 00186 00187 ! create para_env, might need a proper error, bound to this para_env 00188 NULLIFY(para_env_new) 00189 CALL cp_para_env_create(para_env_new,subgroup,error=error) 00190 ! test a sync 00191 CALL mp_sync(para_env_new%group) 00192 00193 ! create blacs, should inherit the preferences for the layout and so on, from the higher level 00194 NULLIFY(blacs_env_new) 00195 CALL cp_blacs_env_create(blacs_env=blacs_env_new, para_env=para_env_new, error=error) 00196 00197 ! create new matrix 00198 NULLIFY(fm_struct_new) 00199 CALL cp_fm_struct_create(fmstruct=fm_struct_new, para_env=para_env_new, context=blacs_env_new,& 00200 nrow_global=n, ncol_global=n, error=error) 00201 CALL cp_fm_create(matrix_new, matrix_struct=fm_struct_new, name="yevd_new_mat", error=error) 00202 CALL cp_fm_create(eigenvectors_new, matrix_struct=fm_struct_new, name="yevd_new_vec", error=error) 00203 00204 ! redistribute old 00205 CALL pdgemr2d(n,n,matrix%local_data(1,1),1,1,matrix%matrix_struct%descriptor, & 00206 matrix_new%local_data(1,1),1,1,matrix_new%matrix_struct%descriptor, & 00207 matrix%matrix_struct%context%group) 00208 00209 ! call scalapack 00210 CALL cp_fm_syevd_base(matrix_new,eigenvectors_new,eig,info,error) 00211 00212 ! redistribute results 00213 CALL pdgemr2d(n,n,eigenvectors_new%local_data(1,1),1,1,eigenvectors_new%matrix_struct%descriptor, & 00214 eigenvectors%local_data(1,1),1,1,eigenvectors%matrix_struct%descriptor, & 00215 eigenvectors%matrix_struct%context%group) 00216 00217 ! free stuff 00218 CALL cp_fm_struct_release(fm_struct_new,error=error) 00219 CALL cp_fm_release(matrix_new,error=error) 00220 CALL cp_fm_release(eigenvectors_new,error=error) 00221 CALL cp_blacs_env_release(blacs_env_new,error=error) 00222 CALL cp_para_env_release(para_env_new,error=error) 00223 00224 ELSE 00225 ! these tasks must help redistribute (they own part of the data), 00226 ! but need fake 'new' data, and their descriptor must indicate this with -1 00227 ! see also scalapack comments on pdgemr2d 00228 fake_descriptor=-1 00229 CALL pdgemr2d(n,n,matrix%local_data(1,1),1,1,matrix%matrix_struct%descriptor, & 00230 fake_local_data(1,1),1,1,fake_descriptor, & 00231 matrix%matrix_struct%context%group) 00232 00233 CALL pdgemr2d(n,n,fake_local_data(1,1),1,1,fake_descriptor, & 00234 eigenvectors%local_data(1,1),1,1,eigenvectors%matrix_struct%descriptor, & 00235 eigenvectors%matrix_struct%context%group) 00236 00237 ! free stuff 00238 CALL mp_comm_free(subgroup) 00239 ENDIF 00240 00241 ! finally, also the eigenvalues need to end up on the non-group member tasks 00242 CALL mp_bcast(eig,0,para_env%group) 00243 00244 ! free more stuff 00245 DEALLOCATE(group_distribution,group_partition) 00246 00247 ELSE 00248 00249 CALL cp_fm_syevd_base(matrix,eigenvectors,eig,info,error) 00250 00251 ENDIF 00252 00253 #else 00254 00255 !MK Retrieve the optimal work array sizes first 00256 myinfo = 0 00257 lwork = -1 00258 liwork = -1 00259 m => matrix%local_data 00260 eig(:) = 0.0_dp 00261 00262 ALLOCATE (work(1),STAT=istat) 00263 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00264 work(:) = 0.0_dp 00265 ALLOCATE (iwork(1),STAT=istat) 00266 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00267 iwork(:) = 0 00268 00269 CALL dsyevd('V','U',n,m(1,1),n,eig(1),work(1),lwork,iwork(1),liwork,myinfo) 00270 00271 IF (myinfo /= 0) THEN 00272 CALL stop_program(routineN,moduleN,__LINE__,& 00273 "ERROR in DSYEVD: Could not retrieve work array sizes") 00274 END IF 00275 00276 ! Reallocate work arrays and perform diagonalisation 00277 lwork = INT(work(1)) 00278 DEALLOCATE (work,STAT=istat) 00279 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00280 ALLOCATE(work(lwork),STAT=istat) 00281 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00282 00283 liwork = iwork(1) 00284 DEALLOCATE (iwork,STAT=istat) 00285 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00286 ALLOCATE(iwork(liwork),STAT=istat) 00287 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00288 iwork(:) = 0 00289 00290 CALL dsyevd('V','U',n,m(1,1),n,eig(1),work(1),lwork,iwork(1),liwork,myinfo) 00291 00292 IF (myinfo /= 0) THEN 00293 CALL stop_program(routineN,moduleN,__LINE__,& 00294 "Matrix diagonalization failed") 00295 END IF 00296 00297 CALL cp_fm_to_fm(matrix,eigenvectors,error=error) 00298 00299 DEALLOCATE (iwork,STAT=istat) 00300 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00301 DEALLOCATE (work,STAT=istat) 00302 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00303 #endif 00304 00305 nmo = SIZE(eigenvalues,1) 00306 IF (nmo > n) THEN 00307 eigenvalues(1:n) = eig(1:n) 00308 ELSE 00309 eigenvalues(1:nmo) = eig(1:nmo) 00310 END IF 00311 00312 DEALLOCATE (eig,STAT=istat) 00313 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00314 CALL timestop(handle) 00315 00316 END SUBROUTINE cp_fm_syevd 00317 00318 SUBROUTINE cp_fm_syevd_base(matrix,eigenvectors,eig,info,error) 00319 00320 TYPE(cp_fm_type), POINTER :: matrix, eigenvectors 00321 REAL(KIND=dp), DIMENSION(:) :: eig 00322 INTEGER, INTENT(OUT), OPTIONAL :: info 00323 TYPE(cp_error_type), INTENT(inout) :: error 00324 00325 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_syevd_base', 00326 routineP = moduleN//':'//routineN 00327 00328 INTEGER :: handle, istat, liwork, lwork, 00329 myinfo, mypcol, myprow, n 00330 INTEGER, DIMENSION(9) :: descm, descv 00331 INTEGER, DIMENSION(:), POINTER :: iwork 00332 LOGICAL :: failure 00333 REAL(KIND=dp), DIMENSION(:), POINTER :: work 00334 REAL(KIND=dp), DIMENSION(:, :), POINTER :: m, v 00335 TYPE(cp_blacs_env_type), POINTER :: context 00336 00337 ! ------------------------------------------------------------------------- 00338 00339 CALL timeset(routineN,handle) 00340 00341 #if defined(__SCALAPACK) 00342 00343 failure = .FALSE. 00344 00345 n = matrix%matrix_struct%nrow_global 00346 m => matrix%local_data 00347 context => matrix%matrix_struct%context 00348 myprow = context%mepos(1) 00349 mypcol = context%mepos(2) 00350 descm(:) = matrix%matrix_struct%descriptor(:) 00351 00352 v => eigenvectors%local_data 00353 descv(:) = eigenvectors%matrix_struct%descriptor(:) 00354 00355 liwork = 7*n + 8*context%num_pe(2) + 2 00356 ALLOCATE(iwork(liwork),STAT=istat) 00357 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00358 ! work space query 00359 00360 lwork = -1 00361 ALLOCATE(work(1),STAT=istat) 00362 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00363 00364 CALL pdsyevd('V','U',n,m(1,1),1,1,descm,eig(1),v(1,1),1,1,descv, & 00365 work(1),lwork,iwork(1),liwork,myinfo) 00366 00367 ! look here for a PDORMTR warning :-) 00368 ! this routine seems to need more workspace than reported 00369 ! (? lapack bug ...) 00370 ! arbitrary additional memory ... we give 100000 more words 00371 ! (it seems to depend on the block size used) 00372 00373 lwork = NINT(work(1)+100000) 00374 ! lwork = NINT(work(1)) 00375 DEALLOCATE (work,STAT=istat) 00376 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00377 ALLOCATE(work(lwork),STAT=istat) 00378 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00379 00380 CALL pdsyevd('V','U',n,m(1,1),1,1,descm,eig(1),v(1,1),1,1,descv, & 00381 work(1),lwork,iwork(1),liwork,myinfo) 00382 00383 IF(.NOT. PRESENT(info) .AND. myinfo /= 0) CALL stop_program & 00384 (routineN,moduleN,__LINE__, "Matrix diagonalization failed") 00385 00386 DEALLOCATE (work,STAT=istat) 00387 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00388 00389 DEALLOCATE (iwork,STAT=istat) 00390 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00391 00392 #endif 00393 00394 IF(PRESENT(info)) info = myinfo 00395 00396 CALL timestop(handle) 00397 00398 END SUBROUTINE cp_fm_syevd_base 00399 00400 ! ***************************************************************************** 00409 SUBROUTINE cp_fm_syevx(matrix,eigenvectors,eigenvalues,neig,work_syevx,error) 00410 00411 ! Diagonalise the symmetric n by n matrix using the LAPACK library. 00412 00413 TYPE(cp_fm_type), POINTER :: matrix 00414 TYPE(cp_fm_type), POINTER, OPTIONAL :: eigenvectors 00415 REAL(KIND = dp), OPTIONAL, INTENT(IN) :: work_syevx 00416 INTEGER, INTENT(IN), OPTIONAL :: neig 00417 REAL(KIND = dp), DIMENSION(:), INTENT(OUT) :: eigenvalues 00418 TYPE(cp_error_type), INTENT(inout) :: error 00419 00420 CHARACTER(LEN=*), PARAMETER :: routineN = "cp_fm_syevx", 00421 routineP = moduleN//":"//routineN 00422 REAL(KIND = dp), PARAMETER :: orfac = -1.0_dp, 00423 vl = 0.0_dp, 00424 vu = 0.0_dp 00425 00426 REAL(KIND = dp) :: abstol, work_syevx_local 00427 INTEGER :: handle,info,istat,liwork,lwork,m,mypcol,myprow,n,nb, 00428 nn,np0,npcol,npe,nprow,nq0,nz,output_unit,itype, neig_local 00429 LOGICAL :: failure, ionode, needs_evecs 00430 00431 INTEGER, DIMENSION(9) :: desca,descz,descb 00432 00433 TYPE(cp_blacs_env_type), POINTER :: context 00434 TYPE(cp_logger_type), POINTER :: logger 00435 REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: gap,w,work 00436 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr,ifail,iwork 00437 REAL(KIND = dp), DIMENSION(:,:), POINTER :: a,z,b 00438 CHARACTER(LEN=1) :: job_type 00439 00440 REAL(KIND = dp), EXTERNAL :: dlamch 00441 00442 #if defined(__SCALAPACK) 00443 INTEGER, EXTERNAL :: iceil,numroc 00444 #else 00445 INTEGER, EXTERNAL :: ilaenv 00446 #endif 00447 ! ------------------------------------------------------------------------- 00448 00449 failure = .FALSE. 00450 ! by default all 00451 n = matrix%matrix_struct%nrow_global 00452 neig_local=n 00453 IF (PRESENT(neig)) neig_local=neig 00454 IF (neig_local == 0) RETURN 00455 00456 CALL timeset(routineN,handle) 00457 00458 needs_evecs=PRESENT(eigenvectors) 00459 00460 logger => cp_error_get_logger(error) 00461 ionode = logger%para_env%mepos==logger%para_env%source 00462 n = matrix%matrix_struct%nrow_global 00463 00464 ! by default allocate all needed space 00465 work_syevx_local=1.0_dp 00466 IF (PRESENT(work_syevx)) work_syevx_local=work_syevx 00467 00468 ! set scalapack job type 00469 IF (needs_evecs) THEN 00470 job_type="V" 00471 ELSE 00472 job_type="N" 00473 ENDIF 00474 00475 ! target the most accurate calculation of the eigenvalues 00476 abstol = 2.0_dp*dlamch("S") 00477 00478 00479 context => matrix%matrix_struct%context 00480 myprow=context%mepos(1) 00481 mypcol=context%mepos(2) 00482 nprow=context%num_pe(1) 00483 npcol=context%num_pe(2) 00484 00485 00486 ALLOCATE (w(n),STAT=istat) 00487 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00488 eigenvalues(:) = 0.0_dp 00489 #if defined(__SCALAPACK) 00490 00491 IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block) THEN 00492 CALL stop_program(routineN,moduleN,__LINE__,& 00493 "Invalid blocksize (no square blocks)") 00494 END IF 00495 00496 a => matrix%local_data 00497 desca(:) = matrix%matrix_struct%descriptor(:) 00498 00499 IF (needs_evecs) THEN 00500 z => eigenvectors%local_data 00501 descz(:) = eigenvectors%matrix_struct%descriptor(:) 00502 ELSE 00503 ! z will not be referenced 00504 z => matrix%local_data 00505 descz=desca 00506 ENDIF 00507 00508 ! Get the optimal work storage size 00509 00510 npe = nprow*npcol 00511 nb = matrix%matrix_struct%nrow_block 00512 nn = MAX(n,nb,2) 00513 np0 = numroc(nn,nb,0,0,nprow) 00514 nq0 = MAX(numroc(nn,nb,0,0,npcol),nb) 00515 00516 IF (needs_evecs) THEN 00517 lwork = 5*n + MAX(5*nn,np0*nq0) + iceil(neig_local,npe)*nn + 2*nb*nb +& 00518 INT(work_syevx_local*REAL((neig_local - 1)*n,dp)) !!!! allocates a full matrix on every CPU !!!!! 00519 ELSE 00520 lwork = 5*n + MAX(5*nn,nb*(np0+1)) 00521 ENDIF 00522 liwork = 6*MAX(N,npe+1,4) 00523 00524 ALLOCATE (gap(npe),STAT=istat) 00525 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00526 gap = 0.0_dp 00527 ALLOCATE (iclustr(2*npe),STAT=istat) 00528 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00529 iclustr = 0 00530 ALLOCATE (ifail(n),STAT=istat) 00531 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00532 ifail = 0 00533 ALLOCATE (iwork(liwork),STAT=istat) 00534 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00535 ALLOCATE (work(lwork),STAT=istat) 00536 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00537 00538 CALL pdsyevx(job_type,"I","U",n,a(1,1),1,1,desca,vl,vu,1,neig_local,abstol,m,nz,w(1),orfac,& 00539 z(1,1),1,1,descz,work(1),lwork,iwork(1),liwork,ifail(1),iclustr(1),gap,info) 00540 00541 ! Error handling 00542 00543 IF (info /= 0) THEN 00544 IF (ionode) THEN 00545 output_unit = cp_logger_get_unit_nr(logger,local=.FALSE.) 00546 WRITE (unit=output_unit,FMT="(/,(T3,A,T12,1X,I10))")& 00547 "info = ",info,& 00548 "lwork = ",lwork,& 00549 "liwork = ",liwork,& 00550 "nz = ",nz 00551 IF (info > 0) THEN 00552 WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,I10)))")& 00553 "ifail = ",ifail 00554 WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,I10)))")& 00555 "iclustr = ",iclustr 00556 WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,E10.3)))")& 00557 "gap = ",gap 00558 END IF 00559 END IF 00560 CALL stop_program(routineN,moduleN,__LINE__,& 00561 "Error in pdsyevx (ScaLAPACK)") 00562 END IF 00563 00564 ! Release work storage 00565 00566 DEALLOCATE (gap,STAT=istat) 00567 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00568 DEALLOCATE (iclustr,STAT=istat) 00569 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00570 DEALLOCATE (ifail,STAT=istat) 00571 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00572 DEALLOCATE (iwork,STAT=istat) 00573 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00574 DEALLOCATE (work,STAT=istat) 00575 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00576 00577 #else 00578 00579 a => matrix%local_data 00580 IF (needs_evecs) THEN 00581 z => eigenvectors%local_data 00582 ELSE 00583 ! z will not be referenced 00584 z => matrix%local_data 00585 ENDIF 00586 00587 ! Get the optimal work storage size 00588 00589 nb = MAX(ilaenv(1,"DSYTRD","U",n,-1,-1,-1),& 00590 ilaenv(1,"DORMTR","U",n,-1,-1,-1)) 00591 00592 lwork = MAX((nb + 3)*n,8*n)+n ! sun bug fix 00593 liwork = 5*n 00594 00595 ALLOCATE (ifail(n),STAT=istat) 00596 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00597 ifail = 0 00598 ALLOCATE (iwork(liwork),STAT=istat) 00599 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00600 ALLOCATE (work(lwork),STAT=istat) 00601 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00602 info = 0 00603 00604 CALL dsyevx(job_type,"I","U",n,a(1,1),n,vl,vu,1,neig_local,abstol,m,w,z(1,1),n,work(1),lwork,& 00605 iwork(1),ifail(1),info) 00606 00607 ! Error handling 00608 00609 IF (info /= 0) THEN 00610 output_unit = cp_logger_get_unit_nr(logger,local=.FALSE.) 00611 WRITE (unit=output_unit,FMT="(/,(T3,A,T12,1X,I10))")& 00612 "info = ",info 00613 IF (info > 0) THEN 00614 WRITE (unit=output_unit,FMT="(/,T3,A,(T12,6(1X,I10)))")& 00615 "ifail = ",ifail 00616 END IF 00617 CALL stop_program(routineN,moduleN,__LINE__,& 00618 "Error in dsyevx") 00619 END IF 00620 00621 ! Release work storage 00622 00623 DEALLOCATE (ifail,STAT=istat) 00624 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00625 DEALLOCATE (iwork,STAT=istat) 00626 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00627 DEALLOCATE (work,STAT=istat) 00628 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00629 00630 #endif 00631 eigenvalues(1:neig_local) = w(1:neig_local) 00632 DEALLOCATE (w,STAT=istat) 00633 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00634 00635 CALL timestop(handle) 00636 00637 END SUBROUTINE cp_fm_syevx 00638 00639 ! ***************************************************************************** 00648 SUBROUTINE cp_fm_syevr(matrix,eigenvectors,eigenvalues,ilow,iup,error) 00649 00650 TYPE(cp_fm_type), POINTER :: matrix 00651 TYPE(cp_fm_type), POINTER, OPTIONAL :: eigenvectors 00652 REAL(KIND = dp), DIMENSION(:), INTENT(OUT) :: eigenvalues 00653 INTEGER, INTENT(IN), OPTIONAL :: ilow,iup 00654 TYPE(cp_error_type), INTENT(inout) :: error 00655 00656 CHARACTER(LEN=*), PARAMETER :: routineN = "cp_fm_syevr", 00657 routineP = moduleN//":"//routineN 00658 REAL(KIND = dp), PARAMETER :: orfac = -1.0_dp, 00659 vl = 0.0_dp, 00660 vu = 0.0_dp 00661 00662 CHARACTER(LEN=1) :: job_type 00663 INTEGER :: handle, ilow_local, info, istat, iup_local, lwork, liwork, m, mypcol,myprow, n, neig, nb, nz 00664 INTEGER, DIMENSION(9) :: desca,descz,descb 00665 INTEGER, DIMENSION(:), ALLOCATABLE :: ifail, iwork 00666 LOGICAL :: failure, ionode, needs_evecs 00667 REAL(dp) :: abstol 00668 REAL(dp), DIMENSION(:), ALLOCATABLE :: w, work 00669 REAL(dp), DIMENSION(:,:), POINTER :: a, z 00670 00671 TYPE(cp_blacs_env_type), POINTER :: context 00672 TYPE(cp_logger_type), POINTER :: logger 00673 00674 REAL(KIND = dp), EXTERNAL :: dlamch 00675 00676 #if defined(__SCALAPACK) 00677 #else 00678 INTEGER, EXTERNAL :: ilaenv 00679 #endif 00680 00681 ! ------------------------------------------------------------------------- 00682 00683 00684 ! by default all 00685 n = matrix%matrix_struct%nrow_global 00686 neig=n 00687 iup_local = n 00688 ilow_local = 1 00689 IF (PRESENT(ilow) .AND. PRESENT(iup)) THEN 00690 neig=iup-ilow+1 00691 iup_local = iup 00692 ilow_local = ilow 00693 END IF 00694 IF (neig <= 0) RETURN 00695 00696 CALL timeset(routineN,handle) 00697 failure =.FALSE. 00698 00699 needs_evecs=PRESENT(eigenvectors) 00700 00701 logger => cp_error_get_logger(error) 00702 ionode = logger%para_env%mepos==logger%para_env%source 00703 n = matrix%matrix_struct%nrow_global 00704 00705 ! set scalapack job type 00706 IF (needs_evecs) THEN 00707 job_type="V" 00708 ELSE 00709 job_type="N" 00710 ENDIF 00711 00712 context => matrix%matrix_struct%context 00713 myprow=context%mepos(1) 00714 mypcol=context%mepos(2) 00715 00716 ALLOCATE(w(n), STAT=istat) 00717 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00718 00719 eigenvalues(:) = 0.0_dp 00720 00721 #if defined(__SCALAPACK) 00722 00723 IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block) THEN 00724 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00725 END IF 00726 00727 a => matrix%local_data 00728 desca(:) = matrix%matrix_struct%descriptor(:) 00729 00730 IF (needs_evecs) THEN 00731 z => eigenvectors%local_data 00732 descz(:) = eigenvectors%matrix_struct%descriptor(:) 00733 ELSE 00734 ! z will not be referenced 00735 z => matrix%local_data 00736 descz=desca 00737 ENDIF 00738 00739 ! First Call: Determine the needed work_space 00740 lwork = -1 00741 ALLOCATE(work(5*n), STAT=istat) 00742 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00743 ALLOCATE(iwork(6*n), STAT=istat) 00744 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00745 #if defined (__SCALAPACK2) 00746 CALL pdsyevr(job_type,'I','U',n,a,1,1,desca,vl,vu,ilow_local,iup_local,m,nz,w(1),z,1,1,descz,work,lwork,iwork,liwork,info) 00747 #endif 00748 lwork = INT(work(1)) 00749 lwork = NINT(work(1)+300000) 00750 liwork = iwork(1) 00751 IF(lwork>SIZE(work,1)) THEN 00752 DEALLOCATE(work,STAT=istat) 00753 ALLOCATE(work(lwork),STAT=istat) 00754 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00755 END IF 00756 IF(liwork>SIZE(iwork,1)) THEN 00757 DEALLOCATE(iwork,STAT=istat) 00758 ALLOCATE(iwork(liwork),STAT=istat) 00759 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00760 END IF 00761 00762 !Second call: solve the eigenvalue problem 00763 info = 0 00764 #if defined (__SCALAPACK2) 00765 CALL pdsyevr(job_type,'I','U',n,a,1,1,desca,vl,vu,ilow_local,iup_local,m,nz,w(1),z,1,1,descz,work,lwork,iwork,liwork,info) 00766 #endif 00767 00768 IF(info>0) THEN 00769 WRITE(*,*) 'Processor ', myprow, mypcol, ': Error! INFO code = ', INFO 00770 END IF 00771 CPPrecondition(info==0,cp_failure_level,routineP,error,failure) 00772 00773 ! Release work storage 00774 DEALLOCATE (iwork,STAT=istat) 00775 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00776 DEALLOCATE (work,STAT=istat) 00777 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00778 00779 #else 00780 00781 a => matrix%local_data 00782 IF (needs_evecs) THEN 00783 z => eigenvectors%local_data 00784 ELSE 00785 ! z will not be referenced 00786 z => matrix%local_data 00787 ENDIF 00788 00789 ! Get the optimal work storage size 00790 00791 nb = MAX(ilaenv(1,"DSYTRD","U",n,-1,-1,-1),& 00792 ilaenv(1,"DORMTR","U",n,-1,-1,-1)) 00793 00794 lwork = MAX((nb + 3)*n,8*n)+n ! sun bug fix 00795 liwork = 5*n 00796 00797 ALLOCATE (ifail(n),STAT=istat) 00798 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00799 ifail = 0 00800 00801 ALLOCATE (iwork(liwork),STAT=istat) 00802 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00803 ALLOCATE (work(lwork),STAT=istat) 00804 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00805 00806 ! target the most accurate calculation of the eigenvalues 00807 abstol = 2.0_dp*dlamch("S") 00808 00809 info = 0 00810 CALL dsyevx(job_type,"I","U",n,a(1,1),n,vl,vu,ilow_local,iup_local,abstol,m,w,z(1,1),n,work(1),lwork,& 00811 iwork(1),ifail(1),info) 00812 00813 ! Error handling 00814 CPPrecondition(info==0,cp_failure_level,routineP,error,failure) 00815 00816 ! Release work storage 00817 DEALLOCATE (iwork,STAT=istat) 00818 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00819 DEALLOCATE (work,STAT=istat) 00820 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00821 00822 #endif 00823 00824 eigenvalues(ilow_local:iup_local) = w(ilow_local:iup_local) 00825 DEALLOCATE (w,STAT=istat) 00826 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00827 00828 CALL timestop(handle) 00829 00830 END SUBROUTINE cp_fm_syevr 00831 00832 ! ***************************************************************************** 00833 SUBROUTINE cp_fm_elpa(matrix,eigenvectors,eigenvalues,error) 00834 00835 TYPE(cp_fm_type), POINTER :: matrix, eigenvectors 00836 REAL(KIND=dp), DIMENSION(:) :: eigenvalues 00837 TYPE(cp_error_type), INTENT(inout) :: error 00838 00839 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_elpa', 00840 routineP = moduleN//':'//routineN 00841 00842 INTEGER :: comm_col, comm_row, group, 00843 handle, istat, mypcol, 00844 myprow, n, n_rows, nblk, neig 00845 INTEGER, DIMENSION(9) :: descm, descv 00846 LOGICAL :: failure 00847 REAL(KIND=dp), DIMENSION(:), POINTER :: eval 00848 REAL(KIND=dp), DIMENSION(:, :), POINTER :: m, m1, v 00849 TYPE(cp_blacs_env_type), POINTER :: context 00850 00851 ! ------------------------------------------------------------------------- 00852 00853 CALL timeset(routineN,handle) 00854 00855 #if defined(__SCALAPACK) 00856 00857 failure = .FALSE. 00858 00859 n = matrix%matrix_struct%nrow_global 00860 m => matrix%local_data 00861 context => matrix%matrix_struct%context 00862 myprow = context%mepos(1) 00863 mypcol = context%mepos(2) 00864 group = matrix%matrix_struct%para_env%group 00865 descm(:) = matrix%matrix_struct%descriptor(:) 00866 00867 00868 ! For ELPA, the MPI communicators along rows/cols are sufficient 00869 ! mpi communicators are created, needed for communicating within 00870 ! rows or columns of processes 00871 CALL mp_comm_split_direct( group, comm_row, mypcol, myprow) 00872 CALL mp_comm_split_direct(group, comm_col, myprow, mypcol) 00873 00874 v => eigenvectors%local_data 00875 descv(:) = eigenvectors%matrix_struct%descriptor(:) 00876 00877 ALLOCATE (m1(SIZE(m,1),SIZE(m,2)),STAT=istat) 00878 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00879 00880 m1 = m 00881 00882 ! elpa needs the full matrix 00883 CALL cp_fm_upper_to_full(matrix,eigenvectors,error=error) 00884 00885 n_rows = matrix%matrix_struct%local_leading_dimension 00886 nblk = matrix%matrix_struct%nrow_block 00887 neig = SIZE(eigenvalues,1) 00888 ! the full eigenvalues vector is needed 00889 ALLOCATE(eval(n),STAT=istat) 00890 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00891 #if defined (__ELPA) 00892 ! Calculate eigenvalues/eigenvectors 00893 CALL solve_evp_real_2stage(n,neig,m,n_rows, eval,v,n_rows,nblk,comm_row,comm_col,group) 00894 #endif 00895 eigenvalues(1:neig) = eval(1:neig) 00896 00897 DEALLOCATE(eval) 00898 DEALLOCATE (m1,STAT=istat) 00899 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00900 00901 ! mpi communicators are freed 00902 CALL mp_comm_free(comm_row) 00903 CALL mp_comm_free(comm_col) 00904 00905 #endif 00906 00907 CALL timestop(handle) 00908 00909 END SUBROUTINE cp_fm_elpa 00910 00911 00912 ! ***************************************************************************** 00913 SUBROUTINE cp_fm_power(matrix,work,exponent,threshold,n_dependent,verbose,error) 00914 00915 ! Raise the real symmetric n by n matrix to the power given by 00916 ! the exponent. All eigenvectors with a corresponding eigenvalue lower 00917 ! than threshold are quenched. result in matrix 00918 00919 ! - Creation (29.03.1999, Matthias Krack) 00920 ! - Parallelised using BLACS and ScaLAPACK (06.06.2001,MK) 00921 00922 TYPE(cp_fm_type), POINTER :: matrix,work 00923 REAL(KIND = dp), INTENT(IN) :: exponent,threshold 00924 INTEGER, INTENT(OUT) :: n_dependent 00925 LOGICAL, INTENT(IN), OPTIONAL :: verbose 00926 TYPE(cp_error_type), INTENT(inout) :: error 00927 00928 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_power', 00929 routineP = moduleN//':'//routineN 00930 00931 INTEGER :: handle,icol_global,icol_local,ipcol,iprow,irow_global, 00932 irow_local,istat,mypcol,myprow,ncol_block,ncol_global,npcol, 00933 nprow,nrow_block,nrow_global 00934 LOGICAL :: failure, my_verbose 00935 REAL(KIND = dp) :: condition_number,f,p 00936 REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: eigenvalues 00937 REAL(KIND = dp), DIMENSION(:,:), POINTER :: eigenvectors 00938 TYPE(cp_blacs_env_type), POINTER :: context 00939 00940 #if defined(__SCALAPACK) 00941 INTEGER, EXTERNAL :: indxg2l,indxg2p 00942 #endif 00943 ! ------------------------------------------------------------------------- 00944 00945 CALL timeset(routineN,handle) 00946 00947 failure = .FALSE. 00948 my_verbose = .FALSE. 00949 IF(PRESENT(verbose)) my_verbose = verbose 00950 00951 context => matrix%matrix_struct%context 00952 myprow=context%mepos(1) 00953 mypcol=context%mepos(2) 00954 nprow=context%num_pe(1) 00955 npcol=context%num_pe(2) 00956 n_dependent = 0 00957 p = 0.5_dp*exponent 00958 00959 nrow_global = matrix%matrix_struct%nrow_global 00960 ncol_global = matrix%matrix_struct%ncol_global 00961 00962 ALLOCATE (eigenvalues(ncol_global),STAT=istat) 00963 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00964 00965 eigenvalues(:) = 0.0_dp 00966 00967 ! Compute the eigenvectors and eigenvalues 00968 00969 #if defined (__ELPA) 00970 CALL cp_fm_elpa(matrix,work,eigenvalues,error=error) 00971 #else 00972 CALL cp_fm_syevd(matrix,work,eigenvalues,error=error) 00973 #endif 00974 00975 #if defined(__SCALAPACK) 00976 nrow_block = work%matrix_struct%nrow_block 00977 ncol_block = work%matrix_struct%ncol_block 00978 00979 eigenvectors => work%local_data 00980 00981 ! Build matrix**exponent with eigenvector quenching 00982 00983 DO icol_global=1,ncol_global 00984 00985 IF (eigenvalues(icol_global) < threshold) THEN 00986 00987 n_dependent = n_dependent + 1 00988 00989 ipcol = indxg2p(icol_global,ncol_block,mypcol,& 00990 work%matrix_struct%first_p_pos(2),npcol) 00991 00992 IF (mypcol == ipcol) THEN 00993 icol_local = indxg2l(icol_global,ncol_block,mypcol,& 00994 work%matrix_struct%first_p_pos(2),npcol) 00995 DO irow_global=1,nrow_global 00996 iprow = indxg2p(irow_global,nrow_block,myprow,& 00997 work%matrix_struct%first_p_pos(1),nprow) 00998 IF (myprow == iprow) THEN 00999 irow_local = indxg2l(irow_global,nrow_block,myprow,& 01000 work%matrix_struct%first_p_pos(1),nprow) 01001 eigenvectors(irow_local,icol_local) = 0.0_dp 01002 END IF 01003 END DO 01004 END IF 01005 01006 ELSE 01007 01008 f = eigenvalues(icol_global)**p 01009 01010 ipcol = indxg2p(icol_global,ncol_block,mypcol,& 01011 work%matrix_struct%first_p_pos(2),npcol) 01012 01013 IF (mypcol == ipcol) THEN 01014 icol_local = indxg2l(icol_global,ncol_block,mypcol,& 01015 work%matrix_struct%first_p_pos(2),npcol) 01016 DO irow_global=1,nrow_global 01017 iprow = indxg2p(irow_global,nrow_block,myprow,& 01018 work%matrix_struct%first_p_pos(1),nprow) 01019 IF (myprow == iprow) THEN 01020 irow_local = indxg2l(irow_global,nrow_block,myprow,& 01021 work%matrix_struct%first_p_pos(1),nprow) 01022 eigenvectors(irow_local,icol_local) =& 01023 f*eigenvectors(irow_local,icol_local) 01024 END IF 01025 END DO 01026 END IF 01027 01028 END IF 01029 01030 END DO 01031 01032 #else 01033 01034 eigenvectors => work%local_data 01035 01036 ! Build matrix**exponent with eigenvector quenching 01037 01038 DO icol_global=1,ncol_global 01039 01040 IF (eigenvalues(icol_global) < threshold) THEN 01041 01042 n_dependent = n_dependent + 1 01043 eigenvectors(1:nrow_global,icol_global) = 0.0_dp 01044 01045 ELSE 01046 01047 f = eigenvalues(icol_global)**p 01048 eigenvectors(1:nrow_global,icol_global) =& 01049 f*eigenvectors(1:nrow_global,icol_global) 01050 01051 END IF 01052 01053 END DO 01054 01055 #endif 01056 CALL cp_fm_syrk("U","N",ncol_global,1.0_dp,work,1,1,0.0_dp,matrix,error=error) 01057 CALL cp_fm_upper_to_full(matrix,work,error=error) 01058 01059 ! Print some warnings/notes 01060 01061 IF (matrix%matrix_struct%para_env%mepos == 0 .AND. my_verbose) THEN 01062 condition_number = ABS(eigenvalues(ncol_global)/eigenvalues(1)) 01063 WRITE (UNIT=cp_logger_get_default_unit_nr(),FMT="(/,(T2,A,ES15.6))")& 01064 "CP_FM_POWER: smallest eigenvalue:",eigenvalues(1),& 01065 "CP_FM_POWER: largest eigenvalue: ",eigenvalues(ncol_global),& 01066 "CP_FM_POWER: condition number: ",condition_number 01067 IF (eigenvalues(1) <= 0.0_dp) THEN 01068 WRITE (UNIT=cp_logger_get_default_unit_nr(),FMT="(/,T2,A)")& 01069 "WARNING: matrix has a negative eigenvalue, tighten EPS_DEFAULT" 01070 END IF 01071 IF (condition_number > 1.0E12_dp) THEN 01072 WRITE (UNIT=cp_logger_get_default_unit_nr(),FMT="(/,T2,A)")& 01073 "WARNING: high condition number => possibly ill-conditioned matrix" 01074 END IF 01075 END IF 01076 01077 DEALLOCATE (eigenvalues,STAT=istat) 01078 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 01079 01080 CALL timestop(handle) 01081 01082 END SUBROUTINE cp_fm_power 01083 01084 ! ***************************************************************************** 01085 SUBROUTINE cp_fm_block_jacobi(matrix,eigenvectors,eigval,thresh,& 01086 start_sec_block) 01087 01088 ! Calculates block diagonalizazion from full symmetric matrix 01089 ! It has its origin in cp_fm_syevx. This routine rotates only elements 01090 ! which are larger than a threshold thresh. 01091 ! start_sec_block is the start of the second block. 01092 ! IT DOES ONLY ONE SWEEP! 01093 01094 ! - Creation (07.10.2002, Martin Fengler) 01095 ! - Cosmetics (05.04.06,MK) 01096 01097 TYPE(cp_fm_type), POINTER :: eigenvectors,matrix 01098 REAL(KIND = dp), DIMENSION(:), INTENT(IN) :: eigval 01099 INTEGER, INTENT(IN) :: start_sec_block 01100 REAL(KIND = dp), INTENT(IN) :: thresh 01101 01102 CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_block_jacobi', 01103 routineP = moduleN//':'//routineN 01104 01105 INTEGER :: handle 01106 REAL(KIND = dp), DIMENSION(:,:), POINTER :: a,ev 01107 01108 REAL(KIND = dp) :: tan_theta,tau,c,s 01109 INTEGER :: I,J,q,p,bs,N,q_loc 01110 REAL(KIND = dp), DIMENSION (:),ALLOCATABLE :: c_ip 01111 01112 #if defined(__SCALAPACK) 01113 TYPE(cp_blacs_env_type), POINTER :: context 01114 01115 INTEGER :: myprow,mypcol,nprow,npcol,ictxt_loc,block_dim_row,block_dim_col, 01116 info,ev_row_block_size,iam,counter,status,ierr 01117 INTEGER :: source,allgrp,mynumrows,ictxt,temp_row,temp_col,temp_nprow, 01118 temp_npcol,mype,npe,a_row,a_col,a_nb,a_mb,iglob,jglob, 01119 i_loc,j_loc 01120 01121 REAL(KIND = dp), DIMENSION(:,:), ALLOCATABLE :: a_loc, ev_loc 01122 REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: diagonale 01123 INTEGER, DIMENSION(9) :: desca,descz,desc_a_block, 01124 desc_ev_loc 01125 01126 INTEGER, EXTERNAL :: numroc,indxl2g,indxg2l 01127 #endif 01128 01129 ! ------------------------------------------------------------------------- 01130 01131 CALL timeset(routineN,handle) 01132 01133 #if defined(__SCALAPACK) 01134 01135 ! Umgebung fuer die Kohn-Sham-Matrix und die Eigenvektoren uebernehmen 01136 01137 context => matrix%matrix_struct%context 01138 source = matrix%matrix_struct%para_env%source 01139 allgrp = matrix%matrix_struct%para_env%group 01140 01141 myprow=context%mepos(1) 01142 mypcol=context%mepos(2) 01143 nprow=context%num_pe(1) 01144 npcol=context%num_pe(2) 01145 01146 N = matrix%matrix_struct%nrow_global 01147 01148 A => matrix%local_data 01149 desca(:) = matrix%matrix_struct%descriptor(:) 01150 EV => eigenvectors%local_data 01151 descz(:) = eigenvectors%matrix_struct%descriptor(:) 01152 01153 ! Kopiere den Block, der wegrotiert werden soll zunaechst auf den Masterprozessor, und anschliessend 01154 ! per Broadcast an alle Prozis 01155 ! ACHTUNG start_sec_block sagt aus WO der ZWEITE Block STARTET!!! 01156 ! Der Block wird mitsamt dem OO-Block bearbeitet 01157 01158 block_dim_row = start_sec_block-1 01159 block_dim_col = N - block_dim_row 01160 ALLOCATE(A_loc(block_dim_row,block_dim_col)) 01161 01162 mype=matrix%matrix_struct%para_env%mepos 01163 npe=matrix%matrix_struct%para_env%num_pe 01164 ictxt =matrix%matrix_struct%context%group 01165 ! get a new context 01166 ictxt_loc=matrix%matrix_struct%para_env%group 01167 CALL cp_blacs_gridinit(ictxt_loc,'Row-major',NPROW*NPCOL,1) 01168 01169 CALL descinit(desc_a_block,block_dim_row,block_dim_col,block_dim_row,& 01170 block_dim_col,0,0,ictxt_loc,block_dim_row,info) 01171 01172 CALL pdgemr2d(block_dim_row,block_dim_col,A,1,start_sec_block,desca,& 01173 A_loc,1,1,desc_a_block,ictxt) 01174 ! Jetzt sind eigentlich nur im Master-Prozess Daten reingekommen 01175 CALL mp_bcast(A_loc,0,allgrp) 01176 01177 ! Da nun jeder ueber den oberen Block verfuegt, koennen wir jetzt die Eigenvektoren so umsortieren, dass 01178 ! man ein NN*1 Prozessgrid hat, und somit jeder Prozessor ueber ein Buendel von Zeilel verfuegt, 01179 ! die selbst unabhaengig modifizieren darf. 01180 01181 ! Aufsetzen der Eigenvektorverteilung 01182 iam = mype 01183 ev_row_block_size = n/(nprow*npcol) 01184 mynumrows = NUMROC(N,ev_row_block_size,iam,0,NPROW*NPCOL) 01185 01186 ALLOCATE(EV_loc(mynumrows,N),c_ip(mynumrows)) 01187 01188 CALL descinit(desc_ev_loc,N,N,ev_row_block_size,N,0,0,ictxt_loc,& 01189 mynumrows,info) 01190 01191 CALL pdgemr2d(N,N,EV,1,1,descz,EV_loc,1,1,desc_ev_loc,ictxt) 01192 01193 ! *** START Diagonalising matrix *** 01194 01195 ! Eigentliche Blockdiagonalisierung 01196 01197 q_loc = 0 01198 01199 DO q=start_sec_block,N 01200 q_loc = q_loc + 1 01201 DO p=1,(start_sec_block-1) 01202 01203 IF(ABS(A_loc(p,q_loc))>thresh) THEN 01204 01205 tau = (eigval(q)-eigval(p))/(2.0_dp*A_loc(p,q_loc)) 01206 01207 tan_theta = SIGN(1.0_dp,tau)/(ABS(tau)+SQRT(1.0_dp+tau*tau)) 01208 01209 ! cos theta 01210 c = 1.0_dp/SQRT(1.0_dp+tan_theta*tan_theta) 01211 s = tan_theta*c 01212 01213 ! Und jetzt noch die Eigenvektoren produzieren: 01214 ! Q * J 01215 ! Verstaendliche Version (bevor die BLAS-Aufrufe sie ersetzt haben) 01216 ! c_ip = c*EV_loc(:,p) - s*EV_loc(:,q) 01217 ! c_iq = s*EV_loc(:,p) + c*EV_loc(:,q) 01218 01219 ! EV(:,p)=c_ip 01220 ! EV(:,q)=c_iq 01221 01222 CALL dcopy(mynumrows,EV_loc(1,p),1,c_ip(1),1) 01223 CALL dscal(mynumrows,c,EV_loc(1,p),1) 01224 CALL daxpy(mynumrows,-s,EV_loc(1,q),1,EV_loc(1,p),1) 01225 CALL dscal(mynumrows,c,EV_loc(1,q),1) 01226 CALL daxpy(mynumrows,s,c_ip(1),1,EV_loc(1,q),1) 01227 01228 END IF 01229 01230 END DO 01231 END DO 01232 01233 ! Nun muessen die Eigenvektoren wieder in die alte Verteilung zurueckverschickt werden. 01234 ! Verschicke EVs zurueck!' 01235 CALL pdgemr2d(N,N,EV_loc,1,1,desc_ev_loc,EV,1,1,descz,ictxt) 01236 01237 ! Speicher freigeben 01238 DEALLOCATE(A_loc,EV_loc,c_ip) 01239 01240 CALL cp_blacs_gridexit(ictxt_loc) 01241 01242 #else 01243 01244 N = matrix%matrix_struct%nrow_global ! Groesse der Matrix A, die bearbeitet werden soll 01245 01246 ALLOCATE(c_ip(N)) ! Speicher fuer den lokalen Eigenwertvektor 01247 01248 A => matrix%local_data ! Contains the Matrix to be worked on 01249 EV => eigenvectors%local_data ! Contains the eigenvectors up to blocksize: Rest ist Muell 01250 01251 ! Start diagonalizing matrix 01252 01253 tan_theta = 0.0_dp 01254 tau = 0.0_dp 01255 01256 DO q=start_sec_block,N 01257 DO p=1,(start_sec_block-1) 01258 01259 IF(ABS(A(p,q))>thresh) THEN 01260 01261 tau = (eigval(q)-eigval(p))/(2.0_dp*A(p,q)) 01262 01263 tan_theta = SIGN(1.0_dp,tau)/(ABS(tau)+SQRT(1.0_dp+tau*tau)) 01264 01265 ! cos theta 01266 c = 1.0_dp/SQRT(1.0_dp+tan_theta*tan_theta) 01267 s = tan_theta*c 01268 01269 ! Und jetzt noch die Eigenvektoren produzieren: 01270 ! Q * J 01271 ! Verstaendliche Version (bevor die BLAS-Aufrufe sie ersetzt haben) 01272 ! c_ip = c*EV(:,p) - s*EV(:,q) 01273 ! c_iq = s*EV(:,p) + c*EV(:,q) 01274 01275 ! EV(:,p)=c_ip 01276 ! EV(:,q)=c_iq 01277 01278 CALL dcopy(N,EV(1,p),1,c_ip(1),1) 01279 CALL dscal(N,c,EV(1,p),1) 01280 CALL daxpy(N,-s,EV(1,q),1,EV(1,p),1) 01281 CALL dscal(N,c,EV(1,q),1) 01282 CALL daxpy(N,s,c_ip(1),1,EV(1,q),1) 01283 01284 END IF 01285 01286 END DO 01287 END DO 01288 01289 ! Release work storage 01290 01291 DEALLOCATE (c_ip) 01292 01293 #endif 01294 01295 CALL timestop(handle) 01296 01297 END SUBROUTINE cp_fm_block_jacobi 01298 01299 END MODULE cp_fm_diag
1.7.3