|
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_dbcsr_operations 00016 USE array_types, ONLY: array_data,& 00017 array_i1d_obj,& 00018 array_new,& 00019 array_nullify,& 00020 array_release 00021 USE cp_blacs_env, ONLY: cp_blacs_env_type,& 00022 get_blacs_info 00023 USE cp_cfm_types, ONLY: cp_cfm_type 00024 USE cp_dbcsr_interface, ONLY: & 00025 cp_dbcsr_add, cp_dbcsr_col_block_sizes, & 00026 cp_dbcsr_complete_redistribute, cp_dbcsr_copy, cp_dbcsr_create, & 00027 cp_dbcsr_distribution, cp_dbcsr_finalize, cp_dbcsr_get_data_size, & 00028 cp_dbcsr_get_data_type, cp_dbcsr_get_info, cp_dbcsr_get_matrix_type, & 00029 cp_dbcsr_get_stored_coordinates, cp_dbcsr_init, & 00030 cp_dbcsr_iterator_blocks_left, cp_dbcsr_iterator_next_block, & 00031 cp_dbcsr_iterator_start, cp_dbcsr_iterator_stop, cp_dbcsr_multiply, & 00032 cp_dbcsr_nblkcols_total, cp_dbcsr_nblkrows_total, cp_dbcsr_norm, & 00033 cp_dbcsr_put_block, cp_dbcsr_release, cp_dbcsr_reserve_block2d, & 00034 cp_dbcsr_reserve_blocks, cp_dbcsr_row_block_sizes, cp_dbcsr_scale, & 00035 cp_dbcsr_valid_index, cp_dbcsr_verify_matrix, cp_dbcsr_work_create 00036 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00037 cp_dbcsr_p_type,& 00038 cp_dbcsr_type 00039 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm 00040 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00041 cp_fm_struct_release,& 00042 cp_fm_struct_type 00043 USE cp_fm_types, ONLY: cp_fm_create,& 00044 cp_fm_get_info,& 00045 cp_fm_release,& 00046 cp_fm_to_fm,& 00047 cp_fm_type 00048 USE cp_para_types, ONLY: cp_blacs_env_type,& 00049 cp_para_env_type 00050 USE dbcsr_block_access, ONLY: dbcsr_access_flush 00051 USE dbcsr_block_operations, ONLY: block_set 00052 USE dbcsr_data_methods, ONLY: dbcsr_get_data,& 00053 dbcsr_get_data_p 00054 USE dbcsr_dist_operations 00055 USE dbcsr_error_handling, ONLY: dbcsr_error_set,& 00056 dbcsr_error_stop,& 00057 dbcsr_error_type 00058 USE dbcsr_io 00059 USE dbcsr_methods, ONLY: & 00060 dbcsr_distribution_col_dist, dbcsr_distribution_init, & 00061 dbcsr_distribution_local_cols, dbcsr_distribution_local_rows, & 00062 dbcsr_distribution_mp, dbcsr_distribution_ncols, & 00063 dbcsr_distribution_new, dbcsr_distribution_nlocal_cols, & 00064 dbcsr_distribution_nlocal_rows, dbcsr_distribution_nrows, & 00065 dbcsr_distribution_release, dbcsr_distribution_row_dist, & 00066 dbcsr_mp_get_coordinates, dbcsr_mp_group, dbcsr_mp_hold, & 00067 dbcsr_mp_mynode, dbcsr_mp_mypcol, dbcsr_mp_myprow, dbcsr_mp_new, & 00068 dbcsr_mp_npcols, dbcsr_mp_nprows, dbcsr_mp_numnodes, dbcsr_mp_pgrid, & 00069 dbcsr_mp_release, dbcsr_wm_use_mutable 00070 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00071 dbcsr_mp_obj,& 00072 dbcsr_norm_frobenius,& 00073 dbcsr_type_antisymmetric,& 00074 dbcsr_type_complex_8,& 00075 dbcsr_type_no_symmetry,& 00076 dbcsr_type_real_8,& 00077 dbcsr_type_symmetric 00078 USE dbcsr_util, ONLY: convert_sizes_to_offsets 00079 USE dbcsr_work_operations, ONLY: add_work_coordinate 00080 USE distribution_2d_types, ONLY: distribution_2d_get,& 00081 distribution_2d_type 00082 USE kinds, ONLY: dp,& 00083 int_size,& 00084 sp 00085 USE message_passing, ONLY: mp_bcast 00086 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 00087 neighbor_list_iterate,& 00088 neighbor_list_iterator_create,& 00089 neighbor_list_iterator_p_type,& 00090 neighbor_list_iterator_release,& 00091 neighbor_list_set_p_type 00092 USE termination, ONLY: stop_memory,& 00093 stop_program 00094 USE timings, ONLY: timeset,& 00095 timestop 00096 00097 !$ USE OMP_LIB 00098 #include "cp_common_uses.h" 00099 00100 IMPLICIT NONE 00101 00102 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_operations' 00103 00104 00105 PUBLIC :: cp_dbcsr_multiply_local, cp_dbcsr_multiply_vec, cp_dbcsr_pack_vec,& 00106 packed_vec_scale, cp_dbcsr_unpack_vec,& 00107 cp_dbcsr_mult_pack_vec_local, packed_vec_bif_tech, & 00108 packed_vec_bif_tech2, cp_dbcsr_copy_vec, packed_vec_build_u, packed_vec_bcast, & 00109 packed_vec_ini 00110 00111 INTERFACE cp_dbcsr_multiply_local 00112 MODULE PROCEDURE cp_dbcsr_multiply_local_d,& 00113 cp_dbcsr_multiply_local_s 00114 END INTERFACE 00115 00116 ! CP2K API emulation 00117 PUBLIC :: cp_dbcsr_add_block_node,& 00118 cp_dbcsr_deallocate_matrix,& 00119 cp_dbcsr_allocate_matrix_set, cp_dbcsr_deallocate_matrix_set,& 00120 cp_dbcsr_from_fm, copy_fm_to_dbcsr, copy_dbcsr_to_fm,& 00121 copy_dbcsr_to_cfm, copy_cfm_to_dbcsr, & 00122 cp_dbcsr_sm_fm_multiply, cp_dbcsr_plus_fm_fm_t,& 00123 cp_dbcsr_get_id_nr,& 00124 cp_dbcsr_alloc_block_from_nbl 00125 00126 ! distribution_2d_type compatibility 00127 PUBLIC :: cp_dbcsr_dist2d_to_dist 00128 00129 PUBLIC :: cp_dbcsr_copy_columns_hack 00130 00131 INTERFACE cp_dbcsr_allocate_matrix_set 00132 MODULE PROCEDURE allocate_dbcsr_matrix_set, allocate_dbcsr_matrix_set_2d 00133 END INTERFACE 00134 00135 INTERFACE cp_dbcsr_deallocate_matrix_set 00136 MODULE PROCEDURE deallocate_dbcsr_matrix_set,& 00137 deallocate_dbcsr_matrix_set_2d 00138 END INTERFACE 00139 00140 INTERFACE cp_dbcsr_plus_fm_fm_t 00141 MODULE PROCEDURE cp_dbcsr_plus_fm_fm_t_native 00142 END INTERFACE 00143 00144 PRIVATE 00145 00146 INTEGER, SAVE, PRIVATE :: last_matrix_id=0 00147 00148 CONTAINS 00149 00150 ! ***************************************************************************** 00154 SUBROUTINE cp_dbcsr_copy_columns_hack(matrix_b, matrix_a,& 00155 ncol, source_start, target_start, para_env, blacs_env, error) 00156 00157 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix_b 00158 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_a 00159 INTEGER, INTENT(IN) :: ncol, source_start, 00160 target_start 00161 TYPE(cp_para_env_type), POINTER :: para_env 00162 TYPE(cp_blacs_env_type), POINTER :: blacs_env 00163 TYPE(cp_error_type), INTENT(INOUT) :: error 00164 00165 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_copy_columns_hack', 00166 routineP = moduleN//':'//routineN 00167 00168 INTEGER :: nfullcols_total, 00169 nfullrows_total 00170 TYPE(cp_fm_struct_type), POINTER :: fm_struct 00171 TYPE(cp_fm_type), POINTER :: fm_matrix_a, fm_matrix_b 00172 00173 NULLIFY(fm_matrix_a, fm_matrix_b, fm_struct) 00174 CALL cp_dbcsr_get_info(matrix_a,nfullrows_total=nfullrows_total,nfullcols_total=nfullcols_total) 00175 CALL cp_fm_struct_create(fm_struct,context=blacs_env,nrow_global=nfullrows_total,& 00176 ncol_global=nfullcols_total,para_env=para_env,error=error) 00177 CALL cp_fm_create(fm_matrix_a,fm_struct,name="fm_matrix_a",error=error) 00178 CALL cp_fm_struct_release(fm_struct,error=error) 00179 00180 CALL cp_dbcsr_get_info(matrix_b,nfullrows_total=nfullrows_total,nfullcols_total=nfullcols_total) 00181 CALL cp_fm_struct_create(fm_struct,context=blacs_env,nrow_global=nfullrows_total,& 00182 ncol_global=nfullcols_total,para_env=para_env,error=error) 00183 CALL cp_fm_create(fm_matrix_b,fm_struct,name="fm_matrix_b",error=error) 00184 CALL cp_fm_struct_release(fm_struct,error=error) 00185 00186 CALL copy_dbcsr_to_fm(matrix_a, fm_matrix_a, error=error) 00187 CALL copy_dbcsr_to_fm(matrix_b, fm_matrix_b, error=error) 00188 00189 CALL cp_fm_to_fm(fm_matrix_a, fm_matrix_b, ncol, source_start, target_start) 00190 00191 CALL copy_fm_to_dbcsr(fm_matrix_b, matrix_b, error=error) 00192 00193 CALL cp_fm_release(fm_matrix_a, error=error) 00194 CALL cp_fm_release(fm_matrix_b, error=error) 00195 00196 END SUBROUTINE cp_dbcsr_copy_columns_hack 00197 00198 00199 ! ***************************************************************************** 00207 SUBROUTINE cp_dbcsr_dist2d_to_dist(dist2d, dist, error, mp_obj) 00208 TYPE(distribution_2d_type), INTENT(IN), 00209 TARGET :: dist2d 00210 TYPE(dbcsr_distribution_obj), 00211 INTENT(OUT) :: dist 00212 TYPE(cp_error_type), INTENT(INOUT) :: error 00213 TYPE(dbcsr_mp_obj), INTENT(IN), OPTIONAL :: mp_obj 00214 00215 INTEGER :: mypcol, myproc, myprow, 00216 numproc 00217 INTEGER, DIMENSION(:), POINTER :: col_dist_data, row_dist_data 00218 INTEGER, DIMENSION(:, :), POINTER :: pgrid 00219 TYPE(array_i1d_obj) :: cd, rd 00220 TYPE(cp_blacs_env_type), POINTER :: blacs_env 00221 TYPE(cp_para_env_type), POINTER :: para_env 00222 TYPE(dbcsr_mp_obj) :: mp_env 00223 TYPE(distribution_2d_type), POINTER :: dist2d_p 00224 00225 ! 00226 ! --------------------------------------------------------------------------- 00227 00228 dist2d_p => dist2d 00229 CALL distribution_2d_get(dist2d_p, error=error,& 00230 row_distribution=row_dist_data, col_distribution=col_dist_data,& 00231 blacs_env=blacs_env) 00232 CALL get_blacs_info(blacs_env, para_env=para_env,& 00233 my_process_row=myprow, my_process_column=mypcol,& 00234 blacs2mpi=pgrid) 00235 myproc = para_env%mepos 00236 numproc = para_env%num_pe 00237 IF (PRESENT (mp_obj)) THEN 00238 mp_env = mp_obj 00239 CALL dbcsr_mp_hold (mp_env) 00240 ELSE 00241 CALL dbcsr_mp_new(mp_env, pgrid, para_env%group, myproc, numproc,& 00242 myprow, mypcol) 00243 ENDIF 00244 CALL array_nullify (rd) 00245 CALL array_nullify (cd) 00246 CALL array_new(rd, row_dist_data) 00247 CALL array_new(cd, col_dist_data) 00248 CALL dbcsr_distribution_new(dist, mp_env, rd, cd) 00249 CALL dbcsr_mp_release (mp_env) 00250 CALL array_release (rd) 00251 CALL array_release (cd) 00252 END SUBROUTINE cp_dbcsr_dist2d_to_dist 00253 00254 ! ***************************************************************************** 00264 00265 SUBROUTINE cp_dbcsr_alloc_block_from_nbl(matrix,sab_orb,error) 00266 00267 TYPE(cp_dbcsr_type) :: matrix 00268 TYPE(neighbor_list_set_p_type), 00269 DIMENSION(:), POINTER :: sab_orb 00270 TYPE(cp_error_type), INTENT(inout) :: error 00271 00272 CHARACTER(LEN=*), PARAMETER :: 00273 routineN = 'cp_dbcsr_alloc_block_from_nbl', 00274 routineP = moduleN//':'//routineN 00275 00276 CHARACTER(LEN=1) :: symmetry 00277 INTEGER :: blk_cnt, handle, iatom, icol, 00278 inode, irow, jatom, last_jatom 00279 INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows, tmp 00280 LOGICAL :: failure, new_atom_b 00281 TYPE(neighbor_list_iterator_p_type), 00282 DIMENSION(:), POINTER :: nl_iterator 00283 00284 CALL timeset(routineN,handle) 00285 00286 failure = .FALSE. 00287 00288 symmetry = cp_dbcsr_get_matrix_type(matrix) 00289 00290 CPPrecondition(ASSOCIATED(sab_orb),cp_failure_level,routineP,error,failure) 00291 00292 CALL cp_dbcsr_finalize (matrix, error=error) 00293 ALLOCATE (rows(1), cols(1)) 00294 blk_cnt = 0 00295 00296 CALL neighbor_list_iterator_create(nl_iterator,sab_orb) 00297 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 00298 CALL get_iterator_info(nl_iterator,iatom=iatom,jatom=jatom,inode=inode) 00299 IF(inode==1) last_jatom = 0 00300 IF (jatom /= last_jatom) THEN 00301 new_atom_b = .TRUE. 00302 last_jatom = jatom 00303 ELSE 00304 new_atom_b = .FALSE. 00305 CYCLE 00306 END IF 00307 IF (blk_cnt+1 .GT. SIZE(rows)) THEN 00308 ALLOCATE (tmp (blk_cnt)) 00309 tmp(:) = rows(:) 00310 DEALLOCATE (rows) 00311 ALLOCATE (rows((blk_cnt+1)*2)) 00312 rows(1:blk_cnt) = tmp(1:blk_cnt) 00313 tmp(:) = cols(:) 00314 DEALLOCATE (cols) 00315 ALLOCATE (cols((blk_cnt+1)*2)) 00316 cols(1:blk_cnt) = tmp(1:blk_cnt) 00317 DEALLOCATE (tmp) 00318 ENDIF 00319 blk_cnt = blk_cnt+1 00320 IF(symmetry==dbcsr_type_no_symmetry) THEN 00321 rows(blk_cnt) = iatom 00322 cols(blk_cnt) = jatom 00323 ELSE 00324 IF(iatom<=jatom) THEN 00325 irow = iatom 00326 icol = jatom 00327 ELSE 00328 irow = jatom 00329 icol = iatom 00330 END IF 00331 rows(blk_cnt) = irow 00332 cols(blk_cnt) = icol 00333 END IF 00334 00335 END DO 00336 CALL neighbor_list_iterator_release(nl_iterator) 00337 00338 ! 00339 CALL cp_dbcsr_reserve_blocks (matrix, rows(1:blk_cnt), cols(1:blk_cnt),& 00340 error=error) 00341 DEALLOCATE (rows) 00342 DEALLOCATE (cols) 00343 CALL cp_dbcsr_finalize( matrix, error=error ) 00344 00345 CALL timestop(handle) 00346 00347 END SUBROUTINE cp_dbcsr_alloc_block_from_nbl 00348 00349 ! ***************************************************************************** 00358 SUBROUTINE cp_dbcsr_multiply_vec(matrix_a, matrix_b, a_row_beg, a_row_end, b_col, pkd, error) 00359 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_a, matrix_b 00360 INTEGER, INTENT(in) :: a_row_beg, a_row_end, b_col 00361 REAL(dp), DIMENSION(:), POINTER :: pkd 00362 TYPE(cp_error_type), INTENT(INOUT) :: error 00363 00364 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_multiply_vec', 00365 routineP = moduleN//':'//routineN 00366 00367 INTEGER :: col, nblkcols_total, node, 00368 pcol, prow, row, timing_handle 00369 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 00370 LOGICAL :: tr 00371 REAL(dp), DIMENSION(:), POINTER :: pkd_b 00372 TYPE(dbcsr_mp_obj) :: mp_obj 00373 00374 CALL timeset(routineN, timing_handle) 00375 00376 nblkcols_total = cp_dbcsr_nblkcols_total(matrix_a) 00377 col_blk_size => array_data (cp_dbcsr_col_block_sizes(matrix_a)) 00378 row_blk_size => array_data (cp_dbcsr_row_block_sizes(matrix_a)) 00379 00380 mp_obj = dbcsr_distribution_mp (cp_dbcsr_distribution (matrix_b)) 00381 00382 00383 ALLOCATE(pkd_b(SIZE(pkd))) 00384 ! 00385 ! packed the b_col 00386 CALL cp_dbcsr_pack_vec(matrix_b, b_col, pkd_b, 'column', error) 00387 ! 00388 ! send the packed col to the right guy 00389 row = 1 00390 col = b_col 00391 !write(*,*) 'cp_dbcsr_multiply_vec: row, col',row,col 00392 tr = .FALSE. 00393 CALL cp_dbcsr_get_stored_coordinates (matrix_b, row, col, tr, node) 00394 !write(*,*) 'cp_dbcsr_multiply_vec: node',node 00395 CALL dbcsr_mp_get_coordinates(mp_obj, node, prow, pcol) 00396 !write(*,*) 'cp_dbcsr_multiply_vec: prow, pcol',prow, pcol 00397 CALL packed_vec_alltoall(pkd_b, pcol, 'column', nblkcols_total, & 00398 row_blk_size*col_blk_size(b_col), mp_obj, error) 00399 ! 00400 ! local multiply 00401 CALL cp_dbcsr_mult_pack_vec_local(matrix_a, pkd_b, a_row_beg, a_row_end, b_col, pkd, error) 00402 ! 00403 ! sum the local products 00404 CALL packed_vec_alltoall(pkd, -1, 'all', nblkcols_total, & 00405 row_blk_size*col_blk_size(b_col), mp_obj, error) 00406 DEALLOCATE(pkd_b) 00407 00408 CALL timestop(timing_handle) 00409 00410 END SUBROUTINE cp_dbcsr_multiply_vec 00411 00412 ! ***************************************************************************** 00421 SUBROUTINE cp_dbcsr_mult_pack_vec_local(matrix_a, pkd_b, a_row_beg, a_row_end, b_col, pkd_c, error) 00422 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_a 00423 REAL(dp), DIMENSION(:), INTENT(in) :: pkd_b 00424 INTEGER, INTENT(in) :: a_row_beg, a_row_end, b_col 00425 REAL(dp), DIMENSION(:), INTENT(inout) :: pkd_c 00426 TYPE(cp_error_type), INTENT(INOUT) :: error 00427 00428 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_mult_pack_vec_local', 00429 routineP = moduleN//':'//routineN 00430 00431 INTEGER :: b_col_size, blk, c_offset, 00432 col, col_size, 00433 nblkrows_total, row, 00434 row_size, timing_handle 00435 INTEGER, DIMENSION(:), POINTER :: col_blk_size 00436 REAL(dp), DIMENSION(:), POINTER :: data_p 00437 TYPE(cp_dbcsr_iterator) :: iter 00438 00439 CALL timeset(routineN, timing_handle) 00440 00441 nblkrows_total = cp_dbcsr_nblkrows_total(matrix_a) 00442 col_blk_size => array_data (cp_dbcsr_col_block_sizes(matrix_a)) 00443 00444 b_col_size = col_blk_size(b_col) 00445 00446 IF (pkd_b(nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 00447 "1pkd_b(n+1).LE.0") 00448 00449 CALL packed_vec_ini(pkd_c, nblkrows_total, error) 00450 !pkd_c(:) = 0.0_dp 00451 c_offset = nblkrows_total + 2 00452 !pkd_c(nblkrows_total+1) = REAL(c_offset,dp) 00453 00454 CALL cp_dbcsr_iterator_start(iter, matrix_a) 00455 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00456 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_p, blk, & 00457 row_size=row_size, col_size=col_size) 00458 IF(row.LT.a_row_beg.OR.row.GT.a_row_end) CYCLE 00459 IF(INT(pkd_b(col)).LE.0) CYCLE 00460 ! 00461 IF(INT(pkd_c(row)).LE.0) THEN 00462 pkd_c(row) = REAL(c_offset,dp) 00463 pkd_c( INT(pkd_c(row)):INT(pkd_c(row))+row_size*b_col_size-1 ) = 0.0_dp 00464 c_offset = c_offset + row_size*b_col_size 00465 pkd_c(nblkrows_total+1) = REAL(c_offset,dp) 00466 ENDIF 00467 CALL dgemm('N','N',row_size,b_col_size,col_size,& 00468 & 1.0_dp,data_p(1),row_size,& 00469 & pkd_b(INT(pkd_b(col))),col_size,& 00470 & 1.0_dp,pkd_c(INT(pkd_c(row))),row_size) 00471 ENDDO 00472 CALL cp_dbcsr_iterator_stop(iter) 00473 00474 IF (pkd_b(nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 00475 "2pkd_b(n+1).LE.0") 00476 IF (pkd_c(nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 00477 "2pkd_c(n+1).LE.0") 00478 00479 CALL timestop(timing_handle) 00480 00481 END SUBROUTINE cp_dbcsr_mult_pack_vec_local 00482 00483 ! ***************************************************************************** 00493 SUBROUTINE cp_dbcsr_multiply_local_d(matrix_a, vec_b, vec_c, ncol, alpha, error) 00494 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_a 00495 REAL(dp), DIMENSION(:, :), INTENT(IN) :: vec_b 00496 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vec_c 00497 INTEGER, INTENT(in), OPTIONAL :: ncol 00498 REAL(dp), INTENT(IN), OPTIONAL :: alpha 00499 TYPE(cp_error_type), INTENT(INOUT) :: error 00500 00501 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_multiply_local_d', 00502 routineP = moduleN//':'//routineN 00503 00504 INTEGER :: blk, col, coloff, my_ncol, 00505 row, rowoff, timing_handle 00506 REAL(dp) :: my_alpha, my_alpha2 00507 REAL(dp), DIMENSION(:, :), POINTER :: data_d 00508 TYPE(cp_dbcsr_iterator) :: iter 00509 00510 CALL timeset(routineN, timing_handle) 00511 00512 00513 my_alpha = 1.0_dp 00514 IF (PRESENT(alpha)) my_alpha = alpha 00515 00516 my_ncol = SIZE(vec_b,2) 00517 IF(PRESENT(ncol)) my_ncol = ncol 00518 00519 my_alpha2 = 0.0_dp 00520 IF(cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_symmetric) my_alpha2 = my_alpha 00521 IF(cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric) my_alpha2 = -my_alpha 00522 00523 CALL cp_dbcsr_iterator_start(iter, matrix_a) 00524 00525 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00526 00527 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_d, blk, row_offset=rowoff, col_offset=coloff) 00528 00529 CALL dgemm('N','N',& 00530 SIZE(data_d,1),my_ncol,SIZE(data_d,2),& 00531 my_alpha, data_d(1,1), SIZE(data_d,1),& 00532 vec_b(coloff,1), SIZE(vec_b,1), & 00533 1.0_dp, vec_c(rowoff,1), SIZE(vec_c,1)) 00534 00535 IF((cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_symmetric.OR.& 00536 cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric)) THEN 00537 IF(row.NE.col) THEN 00538 CALL dgemm('T','N', & 00539 SIZE(data_d,2), my_ncol, SIZE(data_d,1),& 00540 my_alpha2, data_d(1,1), SIZE(data_d,1), & 00541 vec_b(rowoff,1), SIZE(vec_b,1), & 00542 1.0_dp, vec_c(coloff,1), SIZE(vec_c,1)) 00543 ENDIF 00544 ENDIF 00545 ENDDO 00546 00547 CALL cp_dbcsr_iterator_stop(iter) 00548 00549 CALL timestop(timing_handle) 00550 END SUBROUTINE cp_dbcsr_multiply_local_d 00551 00552 SUBROUTINE cp_dbcsr_multiply_local_s(matrix_a, vec_b, vec_c, ncol, alpha, error) 00553 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_a 00554 REAL, DIMENSION(:, :), INTENT(IN) :: vec_b 00555 REAL, DIMENSION(:, :), INTENT(INOUT) :: vec_c 00556 INTEGER, INTENT(in), OPTIONAL :: ncol 00557 REAL, INTENT(IN), OPTIONAL :: alpha 00558 TYPE(cp_error_type), INTENT(INOUT) :: error 00559 00560 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_multiply_local_s', 00561 routineP = moduleN//':'//routineN 00562 00563 INTEGER :: blk, col, coloff, my_ncol, 00564 row, rowoff, timing_handle 00565 REAL :: my_alpha, my_alpha2 00566 REAL, DIMENSION(:, :), POINTER :: data_d 00567 TYPE(cp_dbcsr_iterator) :: iter 00568 00569 CALL timeset(routineN, timing_handle) 00570 00571 my_alpha = 1.0 00572 IF (PRESENT(alpha)) my_alpha = alpha 00573 00574 my_ncol = SIZE(vec_b,2) 00575 IF(PRESENT(ncol)) my_ncol = ncol 00576 00577 my_alpha2 = 0.0 00578 IF(cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_symmetric) my_alpha2 = my_alpha 00579 IF(cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric) my_alpha2 = -my_alpha 00580 00581 CALL cp_dbcsr_iterator_start(iter, matrix_a) 00582 00583 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00584 00585 CALL cp_dbcsr_iterator_next_block(iter, row, col, data_d, blk, row_offset=rowoff, col_offset=coloff) 00586 00587 CALL sgemm('N','N',& 00588 SIZE(data_d,1),my_ncol,SIZE(data_d,2),& 00589 my_alpha, data_d(1,1), SIZE(data_d,1),& 00590 vec_b(coloff,1), SIZE(vec_b,1), & 00591 1.0, vec_c(rowoff,1), SIZE(vec_c,1)) 00592 00593 IF((cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_symmetric.OR.& 00594 cp_dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric)) THEN 00595 IF(row.NE.col) THEN 00596 CALL sgemm('T','N', & 00597 SIZE(data_d,2), my_ncol, SIZE(data_d,1),& 00598 my_alpha2, data_d(1,1), SIZE(data_d,1), & 00599 vec_b(rowoff,1), SIZE(vec_b,1), & 00600 1.0, vec_c(coloff,1), SIZE(vec_c,1)) 00601 ENDIF 00602 ENDIF 00603 ENDDO 00604 00605 CALL cp_dbcsr_iterator_stop(iter) 00606 00607 CALL timestop(timing_handle) 00608 END SUBROUTINE cp_dbcsr_multiply_local_s 00609 00610 00611 00612 ! ***************************************************************************** 00621 SUBROUTINE cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta, error) 00622 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix 00623 TYPE(cp_fm_type), POINTER :: fm_in, fm_out 00624 INTEGER, INTENT(IN) :: ncol 00625 REAL(dp), INTENT(IN), OPTIONAL :: alpha, beta 00626 TYPE(cp_error_type), INTENT(INOUT) :: error 00627 00628 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_sm_fm_multiply', 00629 routineP = moduleN//':'//routineN 00630 00631 INTEGER :: k_in, k_out, timing_handle, 00632 timing_handle_mult 00633 INTEGER, DIMENSION(:), POINTER :: in_col_blk_sizes, 00634 out_col_blk_sizes 00635 TYPE(cp_dbcsr_type) :: in, out 00636 TYPE(array_i1d_obj) :: col_blk_size_right_in, 00637 col_blk_size_right_out 00638 REAL(dp) :: my_alpha, my_beta 00639 TYPE(dbcsr_distribution_obj) :: dist_right_in, product_dist 00640 00641 CALL timeset(routineN, timing_handle) 00642 00643 my_alpha=1.0_dp 00644 my_beta=0.0_dp 00645 IF (PRESENT(alpha)) my_alpha=alpha 00646 IF (PRESENT(beta)) my_beta=beta 00647 00648 CALL cp_fm_get_info(fm_out, ncol_global=k_out, error=error) 00649 00650 CALL cp_fm_get_info(fm_in, ncol_global=k_in, error=error) 00651 !write(*,*)routineN//" -----------------------------------" 00652 !IF (k_in .NE. k_out) & 00653 ! WRITE(*,'(3(A,I5,1X),2(A,F5.2,1X))')& 00654 ! routineN//" ncol", ncol,'k_in',k_in,'k_out',k_out,& 00655 ! 'alpha',my_alpha,'beta',my_beta 00656 00657 IF (ncol.GT.0.AND.k_out.GT.0.AND.k_in.GT.0) THEN 00658 00659 CALL dbcsr_create_dist_r_unrot (dist_right_in, matrix%matrix%m%dist, k_in, & 00660 col_blk_size_right_in) 00661 CALL cp_dbcsr_init(in, error) 00662 CALL cp_dbcsr_create(in, "D", dist_right_in, dbcsr_type_no_symmetry, & 00663 cp_dbcsr_row_block_sizes(matrix), col_blk_size_right_in,& 00664 0, 0, error=error) 00665 00666 CALL cp_dbcsr_init(out, error) 00667 CALL dbcsr_distribution_new (product_dist,& 00668 dbcsr_distribution_mp (cp_dbcsr_distribution(matrix)),& 00669 dbcsr_distribution_row_dist (cp_dbcsr_distribution(matrix)),& 00670 dbcsr_distribution_col_dist (dist_right_in)) 00671 in_col_blk_sizes => array_data (col_blk_size_right_in) 00672 CALL array_nullify (col_blk_size_right_out) 00673 CALL array_new (col_blk_size_right_out, in_col_blk_sizes, lb=1) 00674 out_col_blk_sizes => array_data (col_blk_size_right_out) 00675 CALL match_col_sizes (out_col_blk_sizes, in_col_blk_sizes, k_out) 00676 00677 !if (k_in .ne. k_out) then 00678 ! write(*,*)routineN//" in cs", in_col_blk_sizes 00679 ! write(*,*)routineN//" out cs", out_col_blk_sizes 00680 !endif 00681 00682 CALL cp_dbcsr_create(out, "D", product_dist, dbcsr_type_no_symmetry, & 00683 cp_dbcsr_row_block_sizes(matrix), col_blk_size_right_out,& 00684 0, 0, error=error) 00685 00686 CALL copy_fm_to_dbcsr(fm_in, in, error=error) 00687 IF(ncol.NE.k_out.OR.my_beta.NE.0.0_dp) & 00688 CALL copy_fm_to_dbcsr(fm_out, out, error=error) 00689 00690 CALL timeset(routineN//'_core', timing_handle_mult) 00691 CALL cp_dbcsr_multiply("N", "N", my_alpha, matrix, in,& 00692 my_beta, out, last_column=ncol, error=error) 00693 CALL timestop(timing_handle_mult) 00694 00695 CALL copy_dbcsr_to_fm(out, fm_out,error) 00696 00697 CALL cp_dbcsr_release(in, error=error) 00698 CALL cp_dbcsr_release(out, error=error) 00699 CALL array_release(col_blk_size_right_in) 00700 CALL array_release(col_blk_size_right_out) 00701 CALL dbcsr_distribution_release(dist_right_in) 00702 CALL dbcsr_distribution_release(product_dist) 00703 00704 ENDIF 00705 00706 CALL timestop(timing_handle) 00707 00708 END SUBROUTINE cp_dbcsr_sm_fm_multiply 00709 00710 SUBROUTINE match_col_sizes (sizes1, sizes2, full_num) 00711 INTEGER, DIMENSION(:), INTENT(INOUT) :: sizes1 00712 INTEGER, DIMENSION(:), INTENT(IN) :: sizes2 00713 INTEGER, INTENT(IN) :: full_num 00714 00715 INTEGER :: left, n1, n2, p, rm, used 00716 00717 n1 = SIZE(sizes1) 00718 n2 = SIZE(sizes2) 00719 CALL cp_assert (n1 .EQ. n2, cp_fatal_level, cp_caller_error,& 00720 "match_col_sizes", "distributions must be equal!") 00721 sizes1(1:n1) = sizes2(1:n1) 00722 used = SUM (sizes1(1:n1)) 00723 ! If sizes1 does not cover everything, then we increase the 00724 ! size of the last block; otherwise we reduce the blocks 00725 ! (from the end) until it is small enough. 00726 IF (used .LT. full_num) THEN 00727 sizes1(n1) = sizes1(n1) + full_num-used 00728 ELSE 00729 left = used - full_num 00730 p = n1 00731 DO WHILE (left .GT. 0 .AND. p .GT. 0) 00732 rm = MIN(left, sizes1(p)) 00733 sizes1(p) = sizes1(p) - rm 00734 left = left - rm 00735 p = p-1 00736 ENDDO 00737 ENDIF 00738 END SUBROUTINE match_col_sizes 00739 00740 SUBROUTINE cp_dbcsr_plus_fm_fm_t_native(sparse_matrix,matrix_v,matrix_g,ncol,& 00741 alpha,error) 00742 TYPE(cp_dbcsr_type), INTENT(INOUT) :: sparse_matrix 00743 TYPE(cp_fm_type), POINTER :: matrix_v 00744 TYPE(cp_fm_type), OPTIONAL, POINTER :: matrix_g 00745 INTEGER, INTENT(IN) :: ncol 00746 REAL(KIND=dp), INTENT(IN), OPTIONAL :: alpha 00747 TYPE(cp_error_type), INTENT(inout) :: error 00748 00749 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_plus_fm_fm_t_native', 00750 routineP = moduleN//':'//routineN 00751 00752 INTEGER :: k, nao, timing_handle 00753 LOGICAL :: check_product 00754 REAL(KIND=dp) :: my_alpha, norm 00755 TYPE(array_i1d_obj) :: col_blk_size_left, 00756 col_dist_left 00757 TYPE(cp_dbcsr_type) :: mat_g, mat_v, sparse_matrix2, 00758 sparse_matrix3 00759 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00760 TYPE(cp_fm_type), POINTER :: fm_matrix 00761 TYPE(dbcsr_distribution_obj) :: dist_left 00762 TYPE(dbcsr_mp_obj) :: mp 00763 00764 check_product = .FALSE. 00765 00766 CALL timeset(routineN, timing_handle) 00767 IF (ncol .GT. 0) THEN 00768 CALL cp_assert (cp_dbcsr_valid_index (sparse_matrix), cp_fatal_level,& 00769 cp_caller_error, routineN, "sparse_matrix must pre-exist", error) 00770 ! 00771 ! Setup matrix_v 00772 CALL cp_fm_get_info(matrix_v, ncol_global=k, error=error) 00773 !WRITE(*,*)routineN//'truncated mult k, ncol',k,ncol,' PRESENT (matrix_g)',PRESENT (matrix_g) 00774 mp = dbcsr_distribution_mp (cp_dbcsr_distribution(sparse_matrix)) 00775 CALL create_bl_distribution (col_dist_left, col_blk_size_left,& 00776 k, dbcsr_mp_npcols (mp)) 00777 CALL dbcsr_distribution_new (dist_left, mp,& 00778 dbcsr_distribution_row_dist (cp_dbcsr_distribution(sparse_matrix)),& 00779 col_dist_left) 00780 CALL array_release (col_dist_left) 00781 CALL cp_dbcsr_init (mat_v, error) 00782 CALL cp_dbcsr_create(mat_v, "DBCSR matrix_v", dist_left, dbcsr_type_no_symmetry,& 00783 cp_dbcsr_row_block_sizes (sparse_matrix), col_blk_size_left, 0, 0,& 00784 cp_dbcsr_get_data_type (sparse_matrix), error=error) 00785 CALL copy_fm_to_dbcsr(matrix_v, mat_v, error=error) 00786 CALL cp_dbcsr_verify_matrix(mat_v, error) 00787 ! 00788 ! Setup matrix_g 00789 IF(PRESENT (matrix_g)) THEN 00790 CALL cp_dbcsr_init(mat_g, error) 00791 CALL cp_dbcsr_create(mat_g, "DBCSR matrix_g", dist_left,& 00792 dbcsr_type_no_symmetry,& 00793 cp_dbcsr_row_block_sizes (sparse_matrix),& 00794 cp_dbcsr_col_block_sizes (mat_v),& 00795 data_type=cp_dbcsr_get_data_type (sparse_matrix), error=error) 00796 CALL copy_fm_to_dbcsr(matrix_g, mat_g, error=error) 00797 ENDIF 00798 ! 00799 CALL array_release (col_blk_size_left) 00800 CALL dbcsr_distribution_release (dist_left) 00801 ! 00802 ! 00803 IF(check_product) THEN 00804 NULLIFY(fm_matrix) 00805 CALL cp_fm_get_info(matrix_v,nrow_global=nao,error=error) 00806 CALL cp_fm_struct_create(fm_struct_tmp,context=matrix_v%matrix_struct%context,nrow_global=nao,& 00807 ncol_global=nao,para_env=matrix_v%matrix_struct%para_env,error=error) 00808 CALL cp_fm_create(fm_matrix,fm_struct_tmp,name="fm matrix",error=error) 00809 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00810 CALL copy_dbcsr_to_fm(sparse_matrix,fm_matrix, error=error) 00811 CALL cp_dbcsr_init(sparse_matrix3, error) 00812 CALL cp_dbcsr_copy(sparse_matrix3,sparse_matrix,error=error) 00813 ENDIF 00814 ! 00815 my_alpha = 1.0_dp 00816 IF(PRESENT (alpha)) my_alpha = alpha 00817 IF(PRESENT (matrix_g)) THEN 00818 CALL cp_dbcsr_multiply("N", "T", my_alpha, mat_v, mat_g,& 00819 1.0_dp, sparse_matrix,& 00820 retain_sparsity=.TRUE.,& 00821 last_k = ncol,& 00822 error=error) 00823 ELSE 00824 CALL cp_dbcsr_multiply("N", "T", my_alpha, mat_v, mat_v,& 00825 1.0_dp, sparse_matrix,& 00826 retain_sparsity=.TRUE.,& 00827 last_k = ncol,& 00828 error=error) 00829 ENDIF 00830 00831 IF(check_product) THEN 00832 IF(PRESENT (matrix_g)) THEN 00833 CALL cp_fm_gemm("N","T",nao,nao,ncol,my_alpha,matrix_v,matrix_g,& 00834 1.0_dp,fm_matrix,error=error) 00835 ELSE 00836 CALL cp_fm_gemm("N","T",nao,nao,ncol,my_alpha,matrix_v,matrix_v,& 00837 1.0_dp,fm_matrix,error=error) 00838 ENDIF 00839 00840 CALL cp_dbcsr_init(sparse_matrix2, error) 00841 CALL cp_dbcsr_copy(sparse_matrix2,sparse_matrix,error=error) 00842 CALL cp_dbcsr_scale(sparse_matrix2,alpha_scalar=0.0_dp,error=error) 00843 CALL copy_fm_to_dbcsr(fm_matrix,sparse_matrix2,keep_sparsity=.TRUE., error=error) 00844 CALL cp_dbcsr_add(sparse_matrix2,sparse_matrix,alpha_scalar=1.0_dp,& 00845 beta_scalar=-1.0_dp,error=error) 00846 CALL cp_dbcsr_norm(sparse_matrix2,which_norm=dbcsr_norm_frobenius,& 00847 norm_scalar=norm,error=error) 00848 WRITE(*,*) 'nao=',nao,' k=',k,' ncol=',ncol,' my_alpha=',my_alpha 00849 WRITE(*,*) 'PRESENT (matrix_g)',PRESENT (matrix_g) 00850 WRITE(*,*) 'matrix_type=',cp_dbcsr_get_matrix_type(sparse_matrix) 00851 WRITE(*,*) 'norm(sm+alpha*v*g^t - fm+alpha*v*g^t)/n=',norm/REAL(nao,dp) 00852 IF(norm/REAL(nao,dp).GT.1e-12_dp) THEN 00853 !WRITE(*,*) 'fm_matrix' 00854 !DO j=1,SIZE(fm_matrix%local_data,2) 00855 ! DO i=1,SIZE(fm_matrix%local_data,1) 00856 ! WRITE(*,'(A,I3,A,I3,A,E26.16,A)') 'a(',i,',',j,')=',fm_matrix%local_data(i,j),';' 00857 ! ENDDO 00858 !ENDDO 00859 !WRITE(*,*) 'mat_v' 00860 !CALL cp_dbcsr_print(mat_v,matlab_format=.TRUE.) 00861 !WRITE(*,*) 'mat_g' 00862 !CALL cp_dbcsr_print(mat_g,matlab_format=.TRUE.) 00863 !WRITE(*,*) 'sparse_matrix' 00864 !CALL cp_dbcsr_print(sparse_matrix,matlab_format=.TRUE.) 00865 !WRITE(*,*) 'sparse_matrix2 (-sm + sparse(fm))' 00866 !CALL cp_dbcsr_print(sparse_matrix2,matlab_format=.TRUE.) 00867 !WRITE(*,*) 'sparse_matrix3 (copy of sm input)' 00868 !CALL cp_dbcsr_print(sparse_matrix3,matlab_format=.TRUE.) 00869 !stop 00870 ENDIF 00871 CALL cp_dbcsr_release(sparse_matrix2, error=error) 00872 CALL cp_dbcsr_release(sparse_matrix3, error=error) 00873 CALL cp_fm_release(fm_matrix,error=error) 00874 ENDIF 00875 CALL cp_dbcsr_release (mat_v, error=error) 00876 IF(PRESENT (matrix_g)) CALL cp_dbcsr_release (mat_g, error=error) 00877 ENDIF 00878 CALL timestop(timing_handle) 00879 00880 END SUBROUTINE cp_dbcsr_plus_fm_fm_t_native 00881 00882 ! ***************************************************************************** 00890 SUBROUTINE cp_dbcsr_pack_vec(matrix, ivec, pkd_vec, what, error) 00891 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix 00892 INTEGER, INTENT(IN) :: ivec 00893 REAL(dp), DIMENSION(:), INTENT(INOUT) :: pkd_vec 00894 CHARACTER(LEN=*), INTENT(IN) :: what 00895 TYPE(cp_error_type), INTENT(INOUT) :: error 00896 00897 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_pack_vec', 00898 routineP = moduleN//':'//routineN 00899 00900 INTEGER :: col, col_size, 00901 nblkrows_total, nze, row, 00902 row_size, timing_handle 00903 REAL(KIND=dp), DIMENSION(:), POINTER :: data_p 00904 TYPE(cp_dbcsr_iterator) :: iter 00905 00906 ! --------------------------------------------------------------------------- 00907 00908 CALL timeset(routineN, timing_handle) 00909 00910 nblkrows_total = cp_dbcsr_nblkrows_total(matrix) 00911 pkd_vec(1:nblkrows_total+1) = 0.0_dp ! should be big enough to hold all the data 00912 ! 00913 ! let's go 00914 SELECT CASE(what) 00915 CASE('column') 00916 ! we need nblkrows_total+1 to store rows 00917 ! and 1 extra for adding new data 00918 pkd_vec(nblkrows_total+1) = REAL(nblkrows_total+2,dp) 00919 CALL cp_dbcsr_iterator_start(iter, matrix) 00920 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00921 CALL cp_dbcsr_iterator_next_block (iter, row, col, data_p, & 00922 row_size=row_size, col_size=col_size) 00923 IF(col.ne.ivec) CYCLE 00924 nze = row_size*col_size 00925 !write(*,*) 'row',row,' col',col,' nze',nze 00926 IF(nze.le.0) CYCLE 00927 ! 00928 ! let's copy the block 00929 pkd_vec(row) = pkd_vec(nblkrows_total+1) 00930 pkd_vec(nblkrows_total+1) = pkd_vec(nblkrows_total+1) + REAL(nze,dp) 00931 IF (SIZE(pkd_vec).LT.INT(pkd_vec(nblkrows_total+1))-1) THEN 00932 WRITE(*,*) 'SIZE(pkd_vec)',SIZE(pkd_vec) 00933 WRITE(*,*) 'pkd_vec(nblkrows_total+1)',pkd_vec(nblkrows_total+1) 00934 CALL stop_program(routineN,moduleN,__LINE__,& 00935 "col: SIZE(pkd_vec).LT.pkd_vec(nblkrows_total+1)") 00936 ENDIF 00937 CALL dcopy(nze,data_p(1),1,pkd_vec(INT(pkd_vec(row))),1) 00938 ENDDO 00939 CALL cp_dbcsr_iterator_stop(iter) 00940 CASE('row') 00941 ! we need nblkrows_total+1 to store rows 00942 ! and 1 extra for adding new data 00943 pkd_vec(nblkrows_total+1) = REAL(nblkrows_total+2,dp) 00944 CALL cp_dbcsr_iterator_start(iter, matrix) 00945 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 00946 CALL cp_dbcsr_iterator_next_block (iter, row, col, data_p, & 00947 row_size=row_size, col_size=col_size) 00948 IF(row.ne.ivec) CYCLE 00949 nze = row_size*col_size 00950 !write(*,*) 'row',row,' col',col,' nze',nze 00951 IF(nze.le.0) CYCLE 00952 ! 00953 ! let's copy the block 00954 pkd_vec(col) = pkd_vec(nblkrows_total+1) 00955 pkd_vec(nblkrows_total+1) = pkd_vec(nblkrows_total+1) + nze 00956 IF (SIZE(pkd_vec).LT.INT(pkd_vec(nblkrows_total+1))-1) & 00957 CALL stop_program(routineN,moduleN,__LINE__,& 00958 "SIZE(pkd_vec).LT.pkd_vec(nblkrows_total+1)") 00959 CALL dcopy(nze,data_p(1),1,pkd_vec(INT(pkd_vec(col))),1) 00960 ENDDO 00961 CALL cp_dbcsr_iterator_stop(iter) 00962 CASE DEFAULT 00963 CALL stop_program(routineN,moduleN,__LINE__,"pack what?") 00964 END SELECT 00965 00966 CALL timestop(timing_handle) 00967 00968 END SUBROUTINE cp_dbcsr_pack_vec 00969 00970 ! ***************************************************************************** 00978 SUBROUTINE cp_dbcsr_unpack_vec(matrix, ivec, pkd_vec, what, do_sum, error) 00979 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix 00980 INTEGER, INTENT(IN) :: ivec 00981 REAL(dp), DIMENSION(:), INTENT(IN) :: pkd_vec 00982 CHARACTER(LEN=*), INTENT(IN) :: what 00983 LOGICAL, INTENT(in), OPTIONAL :: do_sum 00984 TYPE(cp_error_type), INTENT(INOUT) :: error 00985 00986 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_unpack_vec', 00987 routineP = moduleN//':'//routineN 00988 00989 INTEGER :: col, col_s, col_size, mynode, nblkcols_total, nblkrows_total, 00990 node, nze, row, row_s, row_size, timing_handle 00991 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 00992 LOGICAL :: my_do_sum, tr 00993 00994 ! --------------------------------------------------------------------------- 00995 00996 CALL timeset(routineN, timing_handle) 00997 00998 my_do_sum = .FALSE. 00999 IF(PRESENT(do_sum))my_do_sum=do_sum 01000 01001 row_blk_size => array_data (cp_dbcsr_row_block_sizes (matrix)) 01002 col_blk_size => array_data (cp_dbcsr_col_block_sizes (matrix)) 01003 nblkrows_total = cp_dbcsr_nblkrows_total(matrix) 01004 nblkcols_total = cp_dbcsr_nblkcols_total(matrix) 01005 mynode = dbcsr_mp_mynode (dbcsr_distribution_mp (cp_dbcsr_distribution (matrix))) 01006 01007 ! 01008 ! let's go 01009 SELECT CASE(what) 01010 CASE('column') 01011 IF (INT(pkd_vec(nblkrows_total+1)).LE.0) & 01012 CALL stop_program(routineN,moduleN,__LINE__,& 01013 "1pkd_vec(n+1).LE.0") 01014 col = ivec 01015 col_size = col_blk_size(col) 01016 DO row = 1, nblkrows_total 01017 IF(INT(pkd_vec(row)).LE.0) CYCLE 01018 row_size = row_blk_size(row) 01019 nze = row_size*col_size 01020 ! add block only on the correct node 01021 tr = .FALSE. 01022 row_s=row ; col_s=col 01023 CALL cp_dbcsr_get_stored_coordinates (matrix, row_s, col_s, tr, node) 01024 IF(node.EQ.mynode) CALL cp_dbcsr_put_block(matrix, row, col, & 01025 pkd_vec(INT(pkd_vec(row)):INT(pkd_vec(row)+nze-1)), & 01026 summation=my_do_sum) 01027 ENDDO ! row 01028 CALL cp_dbcsr_finalize(matrix, error=error) 01029 CASE('row') 01030 IF (INT(pkd_vec(nblkrows_total+1)).LE.0) & 01031 CALL stop_program(routineN,moduleN,__LINE__,& 01032 "2pkd_vec(n+1).LE.0") 01033 row = ivec 01034 row_size = row_blk_size(row) 01035 DO col = 1,nblkcols_total 01036 IF(INT(pkd_vec(col)).LE.0) CYCLE 01037 col_size = col_blk_size(col) 01038 nze = row_size*col_size 01039 ! add block only on the correct node 01040 tr = .FALSE. 01041 row_s=row ; col_s=col 01042 CALL cp_dbcsr_get_stored_coordinates (matrix, row_s, col_s, tr, node) 01043 IF(node.EQ.mynode) CALL cp_dbcsr_put_block(matrix, row, col, & 01044 pkd_vec(INT(pkd_vec(col)):INT(pkd_vec(col)+nze-1)),& 01045 summation=my_do_sum) 01046 ENDDO ! col 01047 CALL cp_dbcsr_finalize(matrix, error=error) 01048 CASE DEFAULT 01049 CALL stop_program(routineN,moduleN,__LINE__,"unpack what?") 01050 END SELECT 01051 01052 CALL timestop(timing_handle) 01053 01054 END SUBROUTINE cp_dbcsr_unpack_vec 01055 01056 ! ***************************************************************************** 01062 SUBROUTINE packed_vec_scale(alpha, block, pkd, nvec, ivec, vec_blk_size, side, error) 01063 REAL(dp), INTENT(IN) :: alpha 01064 REAL(dp), DIMENSION(:), INTENT(in) :: block 01065 REAL(dp), DIMENSION(:), INTENT(inout) :: pkd 01066 INTEGER, INTENT(in) :: nvec, ivec 01067 INTEGER, DIMENSION(:), INTENT(in) :: vec_blk_size 01068 CHARACTER(LEN=*), INTENT(IN) :: side 01069 TYPE(cp_error_type), INTENT(INOUT) :: error 01070 01071 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_scale', 01072 routineP = moduleN//':'//routineN 01073 01074 INTEGER :: i, iblk, k, m, n, offset 01075 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buff 01076 01077 ! --------------------------------------------------------------------------- 01078 01079 ALLOCATE(buff(SIZE(pkd)))!too long buff allocated here... 01080 buff=0.0_dp 01081 SELECT CASE(side) 01082 CASE('left') 01083 iblk = 1 01084 DO i = 1,nvec 01085 offset = INT(pkd(i)) 01086 IF(offset.GT.0) THEN 01087 m = vec_blk_size(i) 01088 n = vec_blk_size(ivec) 01089 k = vec_blk_size(i) 01090 CALL dgemm('N','N',m,n,k,alpha,block(iblk),m,pkd(offset),k,0.0_dp,buff(1),m) 01091 CALL dcopy(m*n,buff(1),1,pkd(offset),1) 01092 ENDIF 01093 iblk = iblk + vec_blk_size(i)**2 01094 ENDDO 01095 CASE('right') 01096 iblk = 1 01097 DO i = 1,nvec 01098 offset = INT(pkd(i)) 01099 IF(offset.GT.0) THEN 01100 m = vec_blk_size(ivec) 01101 n = vec_blk_size(i) 01102 k = vec_blk_size(i) 01103 CALL dgemm('N','N',m,n,k,alpha,pkd(offset),m,block(iblk),k,0.0_dp,buff(1),m) 01104 CALL dcopy(m*n,buff(1),1,pkd(offset),1) 01105 ENDIF 01106 iblk = iblk + vec_blk_size(i)**2 01107 ENDDO 01108 CASE DEFAULT 01109 CALL stop_program(routineN,moduleN,__LINE__,"side?") 01110 END SELECT 01111 DEALLOCATE(buff) 01112 01113 END SUBROUTINE packed_vec_scale 01114 01115 ! ***************************************************************************** 01121 SUBROUTINE packed_vec_bif_tech(mat_v, mat_u, pkd_v_fac, pkd_u_fac, ivec, pkd_v, pkd_u, error) 01122 TYPE(cp_dbcsr_type), INTENT(INOUT) :: mat_v, mat_u 01123 REAL(dp), DIMENSION(:), INTENT(IN) :: pkd_v_fac, pkd_u_fac 01124 INTEGER, INTENT(IN) :: ivec 01125 REAL(dp), DIMENSION(:), INTENT(INOUT) :: pkd_v, pkd_u 01126 TYPE(cp_error_type), INTENT(INOUT) :: error 01127 01128 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_bif_tech', 01129 routineP = moduleN//':'//routineN 01130 01131 INTEGER :: blk, col, col_size, k_row, k_row_size, offset, row, row_size, 01132 timing_handle, u_offset, u_offset_last, ufac_offset, v_offset, 01133 v_offset_last, vfac_offset 01134 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 01135 REAL(KIND=dp), DIMENSION(:), POINTER :: data_u, data_v 01136 01137 ! --------------------------------------------------------------------------- 01138 01139 CALL timeset(routineN, timing_handle) 01140 01141 row_blk_size => array_data (cp_dbcsr_row_block_sizes (mat_u)) 01142 col_blk_size => array_data (cp_dbcsr_col_block_sizes (mat_u)) 01143 01144 IF (pkd_u(cp_dbcsr_nblkrows_total(mat_u)+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01145 "1pkd_u(n+1).LE.0") 01146 IF (pkd_v(cp_dbcsr_nblkrows_total(mat_v)+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01147 "1pkd_v(n+1).LE.0") 01148 CALL dbcsr_get_data (mat_v%matrix, data_v) 01149 CALL dbcsr_get_data (mat_u%matrix, data_u) 01150 ! 01151 !v(k,:) (pkd) = v(k,:) (pkd) - (inv(d(i,i)) * uAk(i,k))' / s (pkd) * v(i,:) (bcsr) i = 1,ivec-1 01152 k_row = ivec 01153 k_row_size = row_blk_size(k_row) 01154 v_offset_last = INT(pkd_v(cp_dbcsr_nblkrows_total(mat_v)+1)) 01155 DO row = 1, cp_dbcsr_nblkrows_total(mat_v) 01156 IF(row.GE.k_row) EXIT !i = 1,ivec-1 01157 row_size = row_blk_size(row) 01158 ufac_offset = INT(pkd_v_fac(row)) 01159 IF(ufac_offset.LE.0) CYCLE 01160 DO blk = mat_v%matrix%m%row_p(row)+1,mat_v%matrix%m%row_p(row+1) 01161 col = mat_v%matrix%m%col_i(blk) 01162 col_size = col_blk_size(col) 01163 offset = ABS(mat_v%matrix%m%blk_p(blk)) 01164 v_offset = INT(pkd_v(col)) 01165 IF(v_offset.LE.0) THEN 01166 v_offset = v_offset_last 01167 v_offset_last = v_offset_last + k_row_size * col_size 01168 pkd_v(v_offset:v_offset+k_row_size * col_size-1) = 0.0_dp 01169 pkd_v(mat_v%matrix%m%nblkrows_total+1) = REAL(v_offset_last,dp) 01170 ENDIF 01171 ! 01172 ! let's multiply and add 01173 CALL dgemm('T', 'N', k_row_size, col_size, row_size,& 01174 & -1.0_dp, pkd_v_fac(ufac_offset), row_size,& 01175 & data_v(offset), row_size,& 01176 & 1.0_dp, pkd_v(v_offset), k_row_size) 01177 ENDDO 01178 ENDDO ! row 01179 ! 01180 !u(k,:) (pkd) = u(k,:) (pkd) - (inv(d(i,i)) * v(i,k))' / s (pkd) * u(i,:) (bcsr) i = 1,ivec-1 01181 u_offset_last = INT(pkd_u(cp_dbcsr_nblkrows_total(mat_u)+1)) 01182 DO row = 1, cp_dbcsr_nblkrows_total(mat_u) 01183 IF(row.GE.k_row) EXIT !i = 1,ivec-1 01184 row_size = row_blk_size(row) 01185 vfac_offset = INT(pkd_u_fac(row)) 01186 IF(vfac_offset.LE.0) CYCLE 01187 DO blk = mat_u%matrix%m%row_p(row)+1,mat_u%matrix%m%row_p(row+1) 01188 col = mat_u%matrix%m%col_i(blk) 01189 col_size = col_blk_size(col) 01190 offset = ABS(mat_u%matrix%m%blk_p(blk)) 01191 u_offset = INT(pkd_u(col)) 01192 IF(u_offset.LE.0) THEN 01193 u_offset = u_offset_last 01194 u_offset_last = u_offset_last + k_row_size * col_size 01195 pkd_u(u_offset:u_offset+k_row_size * col_size-1) = 0.0_dp 01196 pkd_u(mat_u%matrix%m%nblkrows_total+1) = REAL(u_offset_last,dp) 01197 ENDIF 01198 ! 01199 ! let's multiply and add 01200 CALL dgemm('T', 'N', k_row_size, col_size, row_size,& 01201 & -1.0_dp, pkd_u_fac(vfac_offset), row_size,& 01202 & data_u(offset), row_size,& 01203 & 1.0_dp, pkd_u(u_offset), k_row_size) 01204 ENDDO 01205 ENDDO ! row 01206 01207 IF (pkd_u(mat_u%matrix%m%nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01208 "2pkd_u(n+1).LE.0") 01209 IF (pkd_v(mat_v%matrix%m%nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01210 "2pkd_v(n+1).LE.0") 01211 01212 CALL timestop(timing_handle) 01213 01214 END SUBROUTINE packed_vec_bif_tech 01215 01216 01217 SUBROUTINE packed_vec_build_u(pkd_u, pkd_v, k, n, s, vec_blk_size, error) 01218 REAL(dp), DIMENSION(:), INTENT(OUT) :: pkd_u 01219 REAL(dp), DIMENSION(:), INTENT(IN) :: pkd_v 01220 INTEGER, INTENT(IN) :: k, n 01221 REAL(dp), INTENT(IN) :: s 01222 INTEGER, DIMENSION(:), INTENT(IN) :: vec_blk_size 01223 TYPE(cp_error_type), INTENT(INOUT) :: error 01224 01225 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_build_u', 01226 routineP = moduleN//':'//routineN 01227 01228 INTEGER :: i, i_blk_size, k_blk_size, 01229 u_offset 01230 01231 pkd_u(1:n+1) = 0.0_dp 01232 u_offset = n+2 01233 ! 01234 ! u(k,:) = [ v(k,1:k-1)/s -1 0 ... 0] 01235 k_blk_size = vec_blk_size(k) 01236 DO i = 1,k-1 01237 IF(INT(pkd_v(i)).GT.0) THEN 01238 i_blk_size = vec_blk_size(i) 01239 CALL dcopy(k_blk_size*i_blk_size, pkd_v(INT(pkd_v(i))),1,pkd_u(u_offset),1) 01240 CALL dscal(k_blk_size*i_blk_size, -1.0_dp/s, pkd_u(u_offset),1) 01241 pkd_u(i) = u_offset 01242 u_offset = u_offset + k_blk_size*i_blk_size 01243 ENDIF 01244 ENDDO 01245 ! 01246 CALL block_set(k_blk_size, k_blk_size, pkd_u(u_offset:u_offset+ k_blk_size**2-1), & 01247 & 1.0_dp, 0.0_dp) 01248 pkd_u(k) = u_offset 01249 u_offset = u_offset + k_blk_size**2 01250 pkd_u(n+1) = u_offset 01251 01252 END SUBROUTINE packed_vec_build_u 01253 01254 ! ***************************************************************************** 01259 SUBROUTINE packed_vec_bif_tech2(mat_v, pkd_v_fac, ivec, pkd_v, error) 01260 TYPE(cp_dbcsr_type), INTENT(INOUT) :: mat_v 01261 REAL(dp), DIMENSION(:), INTENT(IN) :: pkd_v_fac 01262 INTEGER, INTENT(IN) :: ivec 01263 REAL(dp), DIMENSION(:), INTENT(INOUT) :: pkd_v 01264 TYPE(cp_error_type), INTENT(INOUT) :: error 01265 01266 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_bif_tech2', 01267 routineP = moduleN//':'//routineN 01268 01269 INTEGER :: col, col_size, k_row, 01270 k_row_size, nblkrows_total, 01271 row, row_size, timing_handle 01272 INTEGER, DIMENSION(:), POINTER :: row_blk_size 01273 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buff 01274 REAL(KIND=dp), DIMENSION(:), POINTER :: data_v 01275 TYPE(cp_dbcsr_iterator) :: iter 01276 01277 ! --------------------------------------------------------------------------- 01278 01279 CALL timeset(routineN, timing_handle) 01280 01281 nblkrows_total = cp_dbcsr_nblkrows_total(mat_v) 01282 row_blk_size => array_data (cp_dbcsr_row_block_sizes(mat_v)) 01283 01284 IF (pkd_v(nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01285 "1pkd_v(n+1).LE.0") 01286 01287 ALLOCATE(buff(SIZE(pkd_v))) 01288 01289 CALL packed_vec_ini(buff, nblkrows_total, error) 01290 ! 01291 !v(k,:) (pkd) = v(k,:) (pkd) - (inv(d(i,i)) * uAk(i,k))' / s (pkd) * v(i,:) (bcsr) i = 1,ivec-1 01292 k_row = ivec 01293 k_row_size = row_blk_size(k_row) 01294 01295 CALL cp_dbcsr_iterator_start(iter, mat_v) 01296 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01297 CALL cp_dbcsr_iterator_next_block (iter, row, col, data_v, & 01298 row_size=row_size, col_size=col_size) 01299 IF(row.GE.k_row) EXIT !i = 1,ivec-1 01300 IF(INT(pkd_v_fac(row)).LE.0) CYCLE 01301 IF(INT(pkd_v(col)).LE.0) THEN 01302 pkd_v(col) = pkd_v(nblkrows_total+1) 01303 pkd_v(INT(pkd_v(col)):INT(pkd_v(col)) + k_row_size * col_size-1) = 0.0_dp 01304 !v_offset = v_offset + k_row_size * col_size 01305 pkd_v(nblkrows_total+1) = pkd_v(nblkrows_total+1) + REAL(k_row_size * col_size,dp) 01306 ENDIF 01307 ! 01308 ! let's multiply and add 01309 CALL dgemm('T', 'N', k_row_size, col_size, row_size,& 01310 & -1.0_dp, pkd_v_fac(INT(pkd_v_fac(row))), row_size,& 01311 & data_v(1), row_size,& 01312 & 1.0_dp, pkd_v(INT(pkd_v(col))), k_row_size) 01313 ENDDO 01314 CALL cp_dbcsr_iterator_stop(iter) 01315 ! 01316 IF (pkd_v(nblkrows_total+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01317 "2pkd_v(n+1).LE.0") 01318 DEALLOCATE(buff) 01319 01320 CALL timestop(timing_handle) 01321 01322 END SUBROUTINE packed_vec_bif_tech2 01323 01324 SUBROUTINE packed_vec_ini(pkd_vec, n, error) 01325 REAL(dp), DIMENSION(:), INTENT(OUT) :: pkd_vec 01326 INTEGER, INTENT(IN) :: n 01327 TYPE(cp_error_type), INTENT(INOUT) :: error 01328 01329 pkd_vec(1:n) = 0.0_dp 01330 pkd_vec(n+1) = n + 2 01331 END SUBROUTINE packed_vec_ini 01332 01333 SUBROUTINE packed_vec_alltoall(pkd_vec, source, scope, n, vec_blk_size, mp_obj, error) 01334 REAL(dp), DIMENSION(:), POINTER :: pkd_vec 01335 INTEGER, INTENT(IN) :: source 01336 CHARACTER(LEN=*), INTENT(IN) :: scope 01337 INTEGER, INTENT(IN) :: n 01338 INTEGER, DIMENSION(:), INTENT(in) :: vec_blk_size 01339 TYPE(dbcsr_mp_obj) :: mp_obj 01340 TYPE(cp_error_type), INTENT(INOUT) :: error 01341 01342 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_alltoall', 01343 routineP = moduleN//':'//routineN 01344 01345 INTEGER :: inode, irow, mp_group, 01346 mynode, mypcol, myprow, 01347 npcols, nprows, numnodes, 01348 src, timing_handle 01349 INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi 01350 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buff, buff2 01351 01352 ! --------------------------------------------------------------------------- 01353 01354 CALL timeset(routineN, timing_handle) 01355 01356 numnodes = dbcsr_mp_numnodes (mp_obj) 01357 mynode = dbcsr_mp_mynode (mp_obj) 01358 myprow = dbcsr_mp_myprow (mp_obj) 01359 mypcol = dbcsr_mp_mypcol (mp_obj) 01360 npcols = dbcsr_mp_npcols (mp_obj) 01361 nprows = dbcsr_mp_nprows (mp_obj) 01362 blacs2mpi => dbcsr_mp_pgrid (mp_obj) 01363 mp_group = dbcsr_mp_group (mp_obj) 01364 01365 ALLOCATE(buff(SIZE(pkd_vec))) 01366 ALLOCATE(buff2(10*SIZE(pkd_vec)))! that sucks ! 01367 01368 !buff2(1:int(pkd_vec(n+1))-1) = pkd_vec(1:int(pkd_vec(n+1))-1) 01369 CALL dcopy(INT(pkd_vec(n+1))-1, pkd_vec(1), 1, buff2(1), 1) 01370 01371 SELECT CASE(scope) 01372 CASE('row') 01373 CALL stop_program(routineN,moduleN,__LINE__,"NYI") 01374 CASE('column') 01375 ! 01376 ! simple hack 01377 DO irow = 0,nprows-1 01378 !buff = pkd_vec 01379 CALL dcopy(INT(pkd_vec(n+1))-1, pkd_vec(1), 1, buff(1), 1) 01380 src = blacs2mpi(irow,source) 01381 CALL mp_bcast(buff,src,mp_group) 01382 IF(src.NE.mynode) THEN 01383 CALL packed_vec_add(buff2,buff,vec_blk_size,n,error=error) 01384 ENDIF 01385 ENDDO 01386 CASE('all') 01387 ! 01388 ! simple hack 01389 DO inode = 0,numnodes-1 01390 IF(inode.eq.mynode) THEN 01391 !buff(1:int(pkd_vec(n+1))-1) = pkd_vec(1:int(pkd_vec(n+1))-1) 01392 CALL dcopy(INT(pkd_vec(n+1))-1, pkd_vec(1), 1, buff(1), 1) 01393 !else 01394 !buff(1:int(pkd_vec(n+1))-1) = huge(0.0_dp) 01395 ENDIF 01396 !IF(numnodes.ne.1) then 01397 CALL mp_bcast(buff,inode,mp_group) 01398 IF(inode.ne.mynode) THEN 01399 CALL packed_vec_add(buff2,buff,vec_blk_size,n,error=error) 01400 ENDIF 01401 !endif 01402 ENDDO 01403 CASE DEFAULT 01404 CALL stop_program(routineN,moduleN,__LINE__,"how do you wanna bcast") 01405 END SELECT 01406 01407 IF(SIZE(pkd_vec).lt.INT(buff2(n+1))-1) THEN 01408 DEALLOCATE(pkd_vec) 01409 ALLOCATE(pkd_vec(INT(buff2(n+1)))) 01410 ENDIF 01411 01412 pkd_vec(1:INT(buff2(n+1)) - 1) = buff2(1:INT(buff2(n+1)) - 1) 01413 01414 DEALLOCATE(buff) 01415 DEALLOCATE(buff2) 01416 01417 CALL timestop(timing_handle) 01418 01419 END SUBROUTINE packed_vec_alltoall 01420 01421 SUBROUTINE packed_vec_bcast(pkd_vec, source, scope, do_summation, n, vec_blk_size, mp_obj, error) 01422 REAL(dp), DIMENSION(:), POINTER :: pkd_vec 01423 INTEGER, INTENT(IN) :: source 01424 CHARACTER(LEN=*), INTENT(IN) :: scope 01425 LOGICAL, INTENT(in) :: do_summation 01426 INTEGER, INTENT(IN) :: n 01427 INTEGER, DIMENSION(:), INTENT(in) :: vec_blk_size 01428 TYPE(dbcsr_mp_obj) :: mp_obj 01429 TYPE(cp_error_type), INTENT(INOUT) :: error 01430 01431 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_bcast', 01432 routineP = moduleN//':'//routineN 01433 01434 INTEGER :: icol, inode, irow, mp_group, 01435 mynode, mypcol, myprow, 01436 npcols, nprows, numnodes, 01437 src, timing_handle 01438 INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi 01439 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buff, buff2 01440 01441 ! --------------------------------------------------------------------------- 01442 01443 CALL timeset(routineN, timing_handle) 01444 01445 numnodes = dbcsr_mp_numnodes (mp_obj) 01446 mynode = dbcsr_mp_mynode (mp_obj) 01447 myprow = dbcsr_mp_myprow (mp_obj) 01448 mypcol = dbcsr_mp_mypcol (mp_obj) 01449 npcols = dbcsr_mp_npcols (mp_obj) 01450 nprows = dbcsr_mp_nprows (mp_obj) 01451 blacs2mpi => dbcsr_mp_pgrid (mp_obj) 01452 mp_group = dbcsr_mp_group (mp_obj) 01453 01454 ALLOCATE(buff(SIZE(pkd_vec))) 01455 ALLOCATE(buff2(10*SIZE(pkd_vec)))! that sucks ! 01456 01457 SELECT CASE(scope) 01458 CASE('rowise') 01459 ! 01460 ! simple hack 01461 DO icol = 0,npcols-1 01462 buff = pkd_vec 01463 src = blacs2mpi(source,icol) 01464 CALL mp_bcast(buff,src,mp_group) 01465 IF(mypcol.EQ.icol) THEN 01466 pkd_vec = buff 01467 ENDIF 01468 ENDDO 01469 CASE('columnwise') 01470 ! 01471 ! simple hack 01472 DO irow = 0,nprows-1 01473 buff = pkd_vec 01474 src = blacs2mpi(irow,source) 01475 CALL mp_bcast(buff,src,mp_group) 01476 IF(myprow.EQ.irow) THEN 01477 pkd_vec = buff 01478 ENDIF 01479 ENDDO 01480 CASE('all') 01481 ! 01482 ! simple hack 01483 !buff2(1:int(pkd_vec(n+1))-1) = pkd_vec(1:int(pkd_vec(n+1))-1) 01484 CALL dcopy( INT(pkd_vec(n+1))-1, pkd_vec(1), 1, buff2(1), 1) 01485 DO inode = 0,numnodes-1 01486 IF(inode.eq.mynode) THEN 01487 !buff = pkd_vec 01488 CALL dcopy( INT(pkd_vec(n+1))-1, pkd_vec(1), 1, buff(1), 1) 01489 ELSE 01490 !buff = huge(0.0_dp) 01491 ENDIF 01492 CALL mp_bcast(buff,inode,mp_group) 01493 IF(inode.ne.mynode) THEN 01494 CALL packed_vec_add(buff2,buff,vec_blk_size,n,do_summation,error=error) 01495 ENDIF 01496 ENDDO 01497 IF(SIZE(pkd_vec).lt.INT(buff2(n+1))-1) THEN 01498 DEALLOCATE(pkd_vec) 01499 ALLOCATE(pkd_vec(INT(buff2(n+1)))) 01500 ENDIF 01501 !pkd_vec(1:int(buff2(n+1)) - 1) = buff2(1:int(buff2(n+1)) - 1) 01502 CALL dcopy( INT(buff2(n+1)) - 1, buff2(1), 1, pkd_vec(1), 1) 01503 01504 CASE DEFAULT 01505 CALL stop_program(routineN,moduleN,__LINE__,"how do you wanna bcast") 01506 END SELECT 01507 01508 DEALLOCATE(buff) 01509 DEALLOCATE(buff2) 01510 01511 CALL timestop(timing_handle) 01512 01513 END SUBROUTINE packed_vec_bcast 01514 01515 SUBROUTINE packed_vec_reduce(pkd_vec, to, scope, vec_blk_size, n, mp_obj, error) 01516 REAL(dp), DIMENSION(:), INTENT(INOUT) :: pkd_vec 01517 INTEGER, INTENT(IN) :: to 01518 CHARACTER(LEN=*), INTENT(IN) :: scope 01519 INTEGER, DIMENSION(:), INTENT(IN) :: vec_blk_size 01520 INTEGER, INTENT(IN) :: n 01521 TYPE(dbcsr_mp_obj) :: mp_obj 01522 TYPE(cp_error_type), INTENT(INOUT) :: error 01523 01524 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_reduce', 01525 routineP = moduleN//':'//routineN 01526 01527 INTEGER :: icol, irow, mp_group, mynode, 01528 mypcol, myprow, npcols, 01529 nprows, numnodes, src, 01530 timing_handle 01531 INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi 01532 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buff, sums 01533 01534 ! --------------------------------------------------------------------------- 01535 01536 CALL timeset(routineN, timing_handle) 01537 01538 numnodes = dbcsr_mp_numnodes (mp_obj) 01539 mynode = dbcsr_mp_mynode (mp_obj) 01540 myprow = dbcsr_mp_myprow (mp_obj) 01541 mypcol = dbcsr_mp_mypcol (mp_obj) 01542 npcols = dbcsr_mp_npcols (mp_obj) 01543 nprows = dbcsr_mp_nprows (mp_obj) 01544 blacs2mpi => dbcsr_mp_pgrid (mp_obj) 01545 mp_group = dbcsr_mp_group (mp_obj) 01546 01547 IF (pkd_vec(n+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01548 "a_offset_last.LE.0!!!!!!!") 01549 01550 ALLOCATE(buff(SIZE(pkd_vec)),sums(SIZE(pkd_vec))) 01551 sums = 0.0_dp 01552 SELECT CASE(scope) 01553 CASE('rowise') 01554 CALL stop_program(routineN,moduleN,__LINE__,"more work here") 01555 ! 01556 ! simple hack 01557 DO icol = 0,npcols-1 01558 buff = pkd_vec 01559 src = blacs2mpi(to,icol) 01560 CALL mp_bcast(buff,src,mp_group) 01561 IF(mypcol.EQ.icol) THEN 01562 pkd_vec = buff 01563 CALL packed_vec_add(sums,buff,vec_blk_size,n,error=error) 01564 ENDIF 01565 ENDDO 01566 CASE('columnwise') 01567 ! 01568 ! simple hack 01569 DO irow = 0,nprows-1 01570 buff = pkd_vec 01571 src = blacs2mpi(irow,to) 01572 CALL mp_bcast(buff,src,mp_group) 01573 IF(myprow.EQ.irow) THEN 01574 CALL packed_vec_add(sums,buff,vec_blk_size,n,error=error) 01575 ENDIF 01576 ENDDO 01577 CASE('all') 01578 CALL stop_program(routineN,moduleN,__LINE__,"more work here") 01579 CASE DEFAULT 01580 CALL stop_program(routineN,moduleN,__LINE__,"how do you wanna reduce") 01581 END SELECT 01582 pkd_vec = sums 01583 DEALLOCATE(buff,sums) 01584 01585 IF(pkd_vec(n+1).LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01586 "a_offset_last.LE.0!!!!!!!") 01587 01588 CALL timestop(timing_handle) 01589 01590 END SUBROUTINE packed_vec_reduce 01591 01592 SUBROUTINE packed_vec_add(pkd_vec_a,pkd_vec_b,vec_blk_size,n,do_summation,error) 01593 REAL(dp), DIMENSION(:), INTENT(INOUT) :: pkd_vec_a 01594 REAL(dp), DIMENSION(:), INTENT(IN) :: pkd_vec_b 01595 INTEGER, DIMENSION(:), INTENT(IN) :: vec_blk_size 01596 INTEGER, INTENT(IN) :: n 01597 LOGICAL, INTENT(in), OPTIONAL :: do_summation 01598 TYPE(cp_error_type), INTENT(INOUT) :: error 01599 01600 CHARACTER(LEN=*), PARAMETER :: routineN = 'packed_vec_add', 01601 routineP = moduleN//':'//routineN 01602 01603 INTEGER :: a_offset, a_offset_last, 01604 abeg, aend, b_offset, bbeg, 01605 bend, i, timing_handle 01606 LOGICAL :: my_do_summation 01607 01608 ! --------------------------------------------------------------------------- 01609 ! this points to the last empty entry 01610 01611 CALL timeset(routineN, timing_handle) 01612 01613 my_do_summation = .TRUE. 01614 IF(PRESENT(do_summation)) my_do_summation=do_summation 01615 01616 a_offset_last = INT(pkd_vec_a(n+1)) 01617 IF (a_offset_last.LE.0) CALL stop_program(routineN,moduleN,__LINE__,& 01618 "a_offset_last.LT.0!!!!!!!") 01619 DO i = 1,n 01620 b_offset = INT(pkd_vec_b(i)) 01621 a_offset = INT(pkd_vec_a(i)) 01622 IF(b_offset.LE.0) CYCLE 01623 IF(a_offset.GT.0) THEN 01624 ! the block exsits in a, just add 01625 !CALL daxpy(vec_blk_size(i),1.0_dp,pkd_vec_b(b_offset),1,pkd_vec_a(a_offset),1) 01626 abeg = a_offset 01627 aend = a_offset + vec_blk_size(i) - 1 01628 bbeg = b_offset 01629 bend = b_offset + vec_blk_size(i) - 1 01630 IF(my_do_summation) THEN 01631 pkd_vec_a(abeg:aend) = pkd_vec_a(abeg:aend) + pkd_vec_b(bbeg:bend) 01632 ELSE 01633 pkd_vec_a(abeg:aend) = pkd_vec_b(bbeg:bend) 01634 ENDIF 01635 ELSE 01636 ! the block doesnt exsit in a, copy b at the end 01637 !CALL dcopy(vec_blk_size(i),pkd_vec_b(b_offset),1,pkd_vec_a(a_offset_last),1) 01638 abeg = a_offset_last 01639 aend = a_offset_last + vec_blk_size(i) - 1 01640 bbeg = b_offset 01641 bend = b_offset + vec_blk_size(i) - 1 01642 pkd_vec_a(i) = a_offset_last 01643 pkd_vec_a(abeg:aend) = pkd_vec_b(bbeg:bend) 01644 a_offset_last = a_offset_last + vec_blk_size(i) 01645 ENDIF 01646 ENDDO 01647 ! reset the last empty entry if needed 01648 pkd_vec_a(n+1) = REAL(a_offset_last,dp) 01649 CALL timestop(timing_handle) 01650 END SUBROUTINE packed_vec_add 01651 01652 ! ***************************************************************************** 01660 SUBROUTINE cp_dbcsr_copy_vec(matrix_a, matrix_b, what, ivec, error) 01661 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix_a 01662 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix_b 01663 CHARACTER(LEN=*), INTENT(in) :: what 01664 INTEGER, INTENT(IN) :: ivec 01665 TYPE(cp_error_type), INTENT(INOUT) :: error 01666 01667 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_copy_vec', 01668 routineP = moduleN//':'//routineN 01669 01670 INTEGER :: col, row, timing_handle 01671 LOGICAL :: copy_column 01672 REAL(KIND=dp), DIMENSION(:), POINTER :: data_d 01673 TYPE(cp_dbcsr_iterator) :: iter 01674 01675 ! --------------------------------------------------------------------------- 01676 01677 CALL timeset(routineN, timing_handle) 01678 01679 SELECT CASE(what) 01680 CASE('column') 01681 copy_column = .TRUE. 01682 CASE('row') 01683 copy_column = .FALSE. 01684 CASE DEFAULT 01685 CALL stop_program(routineN,moduleN,__LINE__,"copy what?") 01686 END SELECT 01687 01688 CALL cp_dbcsr_iterator_start(iter, matrix_b) 01689 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01690 CALL cp_dbcsr_iterator_next_block (iter, row, col, data_d) 01691 IF( (copy_column.AND.col.EQ.ivec) .OR. (.NOT.copy_column.AND.row.EQ.ivec) ) & 01692 CALL cp_dbcsr_put_block(matrix_a, row, col, data_d, summation=.FALSE.) 01693 !IF(.NOT.copy_column.AND.row.EQ.ivec) & 01694 ! CALL dbcsr_put_block(matrix_a, row, col, data_d, tr,& 01695 ! summation=.FALSE.) 01696 ENDDO 01697 CALL cp_dbcsr_iterator_stop(iter) 01698 CALL cp_dbcsr_finalize(matrix_a, error=error) 01699 01700 CALL timestop(timing_handle) 01701 01702 END SUBROUTINE cp_dbcsr_copy_vec 01703 01704 01705 ! ***************************************************************************** 01717 SUBROUTINE cp_dbcsr_add_block_node (matrix, block_row, block_col, block, error) 01718 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix 01719 INTEGER, INTENT(IN) :: block_row, block_col 01720 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block 01721 TYPE(cp_error_type), INTENT(INOUT) :: error 01722 01723 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_add_block_node', 01724 routineP = moduleN//':'//routineN 01725 01726 INTEGER :: c, ithread, p, r 01727 LOGICAL :: dbg, existed, is_there, tr 01728 01729 ! --------------------------------------------------------------------------- 01730 ! CALL cp_assert (matrix%m%work_mutable, cp_warning_level, cp_caller_error,& 01731 ! routineN, "Mutable not defined upon DBCSR creation, forcing its use.") 01732 01733 dbg = .FALSE. 01734 01735 ithread = 0 01736 !$ ithread = omp_get_thread_num() 01737 IF (.NOT. ASSOCIATED (matrix%matrix%m%wms)) THEN 01738 CALL cp_dbcsr_work_create (matrix, work_mutable=.TRUE., error=error) 01739 matrix%matrix%m%valid = .FALSE. 01740 ENDIF 01741 !$ CALL cp_assert (SIZE (matrix%matrix%m%wms) .GE. omp_get_num_threads(),& 01742 !$ cp_fatal_level, cp_wrong_args_error, routineN,& 01743 !$ "Too few threads.", error=error) 01744 CALL cp_assert (dbcsr_wm_use_mutable (matrix%matrix%m%wms(ithread+1)),& 01745 cp_warning_level,& 01746 cp_unimplemented_error_nr, routineN,& 01747 "Data loss due to no conversion of appendable to mutable data") 01748 is_there = ASSOCIATED(block) 01749 !r = row ; c = col ; tr = .FALSE. 01750 !CALL dbcsr_get_stored_coordinates (matrix, r, c, tr) 01751 !CALL dbcsr_reserve_block2d (matrix, row, col, block) 01752 !write(*,*) 'add_block_node: block_row',block_row,' block_col',block_col 01753 CALL cp_dbcsr_reserve_block2d (matrix, block_row, block_col, block,& 01754 existed=existed) 01755 ! 01756 IF (dbg) THEN 01757 r = block_row ; c = block_col ; tr = .FALSE. 01758 CALL cp_dbcsr_get_stored_coordinates (matrix, r, c, tr, p) 01759 CALL cp_assert (p .EQ. dbcsr_mp_mynode (dbcsr_distribution_mp (& 01760 cp_dbcsr_distribution(matrix))),& 01761 cp_warning_level, cp_internal_error, routineN,& 01762 "Adding non-local element", error=error) 01763 ENDIF 01764 CALL cp_assert (.NOT.existed, cp_warning_level, cp_wrong_args_error,& 01765 routineN, "You should not add existing blocks according to old API.") 01766 IF(.NOT.is_there) block(:,:) = 0.0_dp 01767 END SUBROUTINE cp_dbcsr_add_block_node 01768 01769 ! ***************************************************************************** 01775 FUNCTION cp_dbcsr_get_id_nr(matrix,error) RESULT(res) 01776 TYPE(cp_dbcsr_type), POINTER :: matrix 01777 TYPE(cp_error_type), INTENT(inout) :: error 01778 INTEGER :: res 01779 01780 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_get_id_nr', 01781 routineP = moduleN//':'//routineN 01782 01783 LOGICAL :: failure 01784 01785 failure=.FALSE. 01786 ! CPPrecondition(associated(matrix),cp_failure_level,routineP,error,failure) 01787 CPPrecondition(matrix%matrix%m%refcount>0,cp_failure_level,routineP,error,failure) 01788 IF (.NOT. failure) THEN 01789 res=matrix%matrix%m%id_nr 01790 ELSE 01791 res=0 01792 END IF 01793 END FUNCTION cp_dbcsr_get_id_nr 01794 01795 ! ***************************************************************************** 01800 SUBROUTINE cp_dbcsr_deallocate_matrix(matrix, error) 01801 TYPE(cp_dbcsr_type), POINTER :: matrix 01802 TYPE(cp_error_type), INTENT(INOUT) :: error 01803 01804 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_deallocate_matrix', 01805 routineP = moduleN//':'//routineN 01806 01807 INTEGER :: timing_handle 01808 01809 CALL timeset(routineN, timing_handle) 01810 01811 CALL cp_dbcsr_release (matrix, error=error) 01812 CALL cp_assert (.NOT. cp_dbcsr_valid_index(matrix), cp_warning_level,& 01813 cp_caller_error, routineN,& 01814 'You should not "deallocate" a referenced matrix. '//& 01815 'Avoid pointers to DBCSR matrices.') 01816 DEALLOCATE (matrix) 01817 NULLIFY (matrix) 01818 01819 CALL timestop(timing_handle) 01820 01821 END SUBROUTINE cp_dbcsr_deallocate_matrix 01822 01823 ! ***************************************************************************** 01824 ! CP2k-compatible matrix sets 01825 ! ***************************************************************************** 01826 01827 01828 ! ***************************************************************************** 01836 SUBROUTINE allocate_dbcsr_matrix_set(matrix_set, nmatrix, error) 01837 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01838 POINTER :: matrix_set 01839 INTEGER, INTENT(IN) :: nmatrix 01840 TYPE(cp_error_type), INTENT(INOUT) :: error 01841 01842 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_dbcsr_matrix_set', 01843 routineP = moduleN//':'//routineN 01844 01845 INTEGER :: imatrix, istat 01846 01847 IF (ASSOCIATED(matrix_set)) CALL cp_dbcsr_deallocate_matrix_set(matrix_set,error=error) 01848 ALLOCATE (matrix_set(nmatrix),STAT=istat) 01849 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01850 "matrix_set",int_size*nmatrix) 01851 DO imatrix=1,nmatrix 01852 NULLIFY (matrix_set(imatrix)%matrix) 01853 !ALLOCATE (matrix_set(imatrix)%matrix) 01854 !CALL cp_dbcsr_init (matrix_set(imatrix)%matrix) 01855 END DO 01856 END SUBROUTINE allocate_dbcsr_matrix_set 01857 01858 ! ***************************************************************************** 01866 SUBROUTINE allocate_dbcsr_matrix_set_2d(matrix_set,nmatrix,mmatrix,error) 01867 TYPE(cp_dbcsr_p_type), DIMENSION(:, :), 01868 POINTER :: matrix_set 01869 INTEGER, INTENT(IN) :: nmatrix, mmatrix 01870 TYPE(cp_error_type), INTENT(inout) :: error 01871 01872 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_dbcsr_matrix_set_2d', 01873 routineP = moduleN//':'//routineN 01874 01875 INTEGER :: imatrix, istat, jmatrix 01876 01877 IF (ASSOCIATED(matrix_set)) CALL cp_dbcsr_deallocate_matrix_set(matrix_set,error=error) 01878 ALLOCATE (matrix_set(nmatrix,mmatrix),STAT=istat) 01879 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01880 "matrix_set",int_size*nmatrix*mmatrix) 01881 DO jmatrix=1,mmatrix 01882 DO imatrix=1,nmatrix 01883 NULLIFY (matrix_set(imatrix,jmatrix)%matrix) 01884 !ALLOCATE (matrix_set(imatrix,jmatrix)%matrix) 01885 !CALL cp_dbcsr_init (matrix_set(imatrix,jmatrix)%matrix) 01886 END DO 01887 END DO 01888 END SUBROUTINE allocate_dbcsr_matrix_set_2d 01889 01890 01891 ! ***************************************************************************** 01898 SUBROUTINE deallocate_dbcsr_matrix_set(matrix_set,error) 01899 01900 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01901 POINTER :: matrix_set 01902 TYPE(cp_error_type), INTENT(inout) :: error 01903 01904 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_dbcsr_matrix_set', 01905 routineP = moduleN//':'//routineN 01906 01907 INTEGER :: imatrix, istat, timing_handle 01908 01909 CALL timeset(routineN, timing_handle) 01910 01911 IF (ASSOCIATED(matrix_set)) THEN 01912 DO imatrix=1,SIZE(matrix_set) 01913 CALL cp_dbcsr_deallocate_matrix(matrix_set(imatrix)%matrix,error=error) 01914 END DO 01915 DEALLOCATE (matrix_set,STAT=istat) 01916 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"matrix_set") 01917 END IF 01918 01919 CALL timestop(timing_handle) 01920 01921 END SUBROUTINE deallocate_dbcsr_matrix_set 01922 01923 ! ***************************************************************************** 01930 SUBROUTINE deallocate_dbcsr_matrix_set_2d(matrix_set,error) 01931 01932 TYPE(cp_dbcsr_p_type), DIMENSION(:, :), 01933 POINTER :: matrix_set 01934 TYPE(cp_error_type), INTENT(inout) :: error 01935 01936 CHARACTER(LEN=*), PARAMETER :: 01937 routineN = 'deallocate_dbcsr_matrix_set_2d', 01938 routineP = moduleN//':'//routineN 01939 01940 INTEGER :: imatrix, istat, jmatrix 01941 01942 IF (ASSOCIATED(matrix_set)) THEN 01943 DO jmatrix=1,SIZE(matrix_set,2) 01944 DO imatrix=1,SIZE(matrix_set,1) 01945 CALL cp_dbcsr_deallocate_matrix(matrix_set(imatrix,jmatrix)%matrix,error=error) 01946 END DO 01947 END DO 01948 DEALLOCATE (matrix_set,STAT=istat) 01949 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"matrix_set") 01950 END IF 01951 END SUBROUTINE deallocate_dbcsr_matrix_set_2d 01952 01953 ! ***************************************************************************** 01958 SUBROUTINE cp_dbcsr_type_from_sm_symmetry (dbcsr_matrix_type, sm_symmetry) 01959 CHARACTER, INTENT(OUT) :: dbcsr_matrix_type 01960 CHARACTER(LEN=*), INTENT(IN) :: sm_symmetry 01961 01962 CHARACTER(LEN=*), PARAMETER :: 01963 routineN = 'cp_dbcsr_type_from_sm_symmetry', 01964 routineP = moduleN//':'//routineN 01965 01966 ! --------------------------------------------------------------------------- 01967 01968 SELECT CASE (sm_symmetry) 01969 CASE ("symmetric") 01970 dbcsr_matrix_type = dbcsr_type_symmetric 01971 CASE ("antisymmetric") 01972 dbcsr_matrix_type = dbcsr_type_antisymmetric 01973 CASE ("none", "no symmetry") 01974 dbcsr_matrix_type = dbcsr_type_no_symmetry 01975 CASE default 01976 CALL cp_assert (.FALSE., cp_warning_level, cp_caller_error, routineN,& 01977 "Unknown matrix type "//sm_symmetry) 01978 END SELECT 01979 END SUBROUTINE cp_dbcsr_type_from_sm_symmetry 01980 01981 01982 ! ***************************************************************************** 01992 SUBROUTINE cp_dbcsr_from_fm(matrix, fm, threshold, distribution, row_blk_size,& 01993 col_blk_size, error) 01994 TYPE(cp_dbcsr_type), INTENT(OUT) :: matrix 01995 TYPE(cp_fm_type), POINTER :: fm 01996 REAL(KIND=dp), INTENT(IN) :: threshold 01997 TYPE(dbcsr_distribution_obj) :: distribution 01998 TYPE(array_i1d_obj), INTENT(IN) :: row_blk_size, col_blk_size 01999 TYPE(cp_error_type), INTENT(INOUT) :: error 02000 02001 CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_dbcsr_from_fm', 02002 routineP = moduleN//':'//routineN 02003 02004 INTEGER :: timing_handle 02005 02006 CALL timeset(routineN, timing_handle) 02007 02008 !CALL cp_dbcsr_init (matrix, error)! the matrix should already be initialized 02009 CALL cp_dbcsr_create(matrix, fm%name, distribution, dbcsr_type_no_symmetry,& 02010 row_blk_size, col_blk_size,& 02011 0, 0, dbcsr_type_real_8, error=error) 02012 CALL copy_fm_to_dbcsr(fm, matrix, error=error) 02013 CALL cp_dbcsr_verify_matrix(matrix, error) 02014 CALL timestop(timing_handle) 02015 02016 END SUBROUTINE cp_dbcsr_from_fm 02017 02018 02019 ! ***************************************************************************** 02036 SUBROUTINE copy_fm_to_dbcsr(fm,matrix,alpha,beta,keep_sparsity,error) 02037 TYPE(cp_fm_type), POINTER :: fm 02038 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix 02039 REAL(kind=dp), INTENT(IN), OPTIONAL :: alpha, beta 02040 LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity 02041 TYPE(cp_error_type), INTENT(inout) :: error 02042 02043 CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_fm_to_dbcsr', 02044 routineP = moduleN//':'//routineN 02045 02046 INTEGER :: blk_p, col, col_l, col_size, error_handler, group, handle, 02047 nblkcols_local, nblkcols_total, nblkrows_local, nblkrows_total, 02048 ncol_block, ncol_global, nfullcols_local, nfullcols_total, 02049 nfullrows_local, nfullrows_total, nrow_block, nrow_global, nze, row, 02050 row_l, row_size 02051 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, 02052 last_col, last_row, 02053 local_col_sizes, 02054 local_row_sizes 02055 INTEGER, DIMENSION(:), POINTER :: cbs, local_cols, local_rows, 02056 rbs 02057 REAL(kind=dp) :: my_beta 02058 REAL(KIND=dp), DIMENSION(:), POINTER :: blk_1d_dp 02059 REAL(KIND=dp), DIMENSION(:, :), POINTER :: fm_block 02060 REAL(KIND=sp), DIMENSION(:), POINTER :: blk_1d_sp 02061 REAL(kind=sp), DIMENSION(:, :), POINTER :: fm_block_sp 02062 TYPE(array_i1d_obj) :: col_blk_size, row_blk_size 02063 TYPE(cp_blacs_env_type), POINTER :: context 02064 TYPE(cp_dbcsr_type) :: bc_mat 02065 TYPE(dbcsr_distribution_obj) :: bc_dist 02066 TYPE(dbcsr_error_type) :: dbcsr_error 02067 02068 ! --------------------------------------------------------------------------- 02069 02070 CALL dbcsr_error_set(routineN, error_handler, dbcsr_error) 02071 CALL dbcsr_access_flush (matrix%matrix, error=dbcsr_error) 02072 CALL timeset(routineN,handle) 02073 02074 my_beta=0._dp 02075 IF (PRESENT(beta)) THEN 02076 CALL cp_assert (beta .EQ. my_beta, cp_fatal_level,& 02077 cp_unimplemented_error_nr, routineN,& 02078 "beta not supported, use matrix addition instead") 02079 my_beta=beta 02080 ENDIF 02081 group = fm%matrix_struct%para_env%group 02082 context => fm%matrix_struct%context 02083 nrow_block = fm%matrix_struct%nrow_block 02084 ncol_block = fm%matrix_struct%ncol_block 02085 nrow_global = fm%matrix_struct%nrow_global 02086 ncol_global = fm%matrix_struct%ncol_global 02087 02088 CALL cp_dbcsr_get_info(matrix,& 02089 nfullrows_total=nfullrows_total,& 02090 nfullcols_total=nfullcols_total) 02091 02092 CALL cp_assert (nrow_global.eq.nfullrows_total,& 02093 cp_fatal_level, cp_caller_error, routineN,& 02094 "FM and DBCSR matrix sizes do not match in rows") 02095 CALL cp_assert (ncol_global.eq.nfullcols_total,& 02096 cp_fatal_level, cp_caller_error, routineN,& 02097 "FM and DBCSR matrix sizes do not match in columns") 02098 02099 ! Create a block-cyclic distribution compatible with the FM matrix. 02100 CALL dbcsr_distribution_init (bc_dist) 02101 CALL dbcsr_create_dist_block_cyclic (bc_dist,& 02102 nfullrows_total, nfullcols_total,& ! Actual full matrix size 02103 nrow_block, ncol_block,& ! BLACS parameters 02104 dbcsr_distribution_mp (cp_dbcsr_distribution (matrix)),& 02105 row_blk_size, col_blk_size) ! block-cyclic row/col sizes 02106 02107 ! Create the block-cyclic DBCSR matrix 02108 CALL cp_dbcsr_init (bc_mat, error) 02109 CALL cp_dbcsr_create (bc_mat, "Block-cyclic "//matrix%matrix%m%name, bc_dist,& 02110 cp_dbcsr_get_matrix_type(matrix), row_blk_size, col_blk_size, 0, 0,& 02111 data_type=cp_dbcsr_get_data_type(matrix),error=error) 02112 02113 !call dbcsr_finalize (bc_mat) 02114 CALL dbcsr_distribution_release (bc_dist) 02115 CALL array_release (row_blk_size) 02116 CALL array_release (col_blk_size) 02117 02118 CALL cp_dbcsr_get_info(bc_mat,& 02119 nblkrows_total=nblkrows_total,& 02120 nblkcols_total=nblkcols_total,& 02121 nblkrows_local=nblkrows_local,& 02122 nblkcols_local=nblkcols_local,& 02123 nfullrows_local=nfullrows_local,& 02124 nfullcols_local=nfullcols_local,& 02125 nfullrows_total=nfullrows_total,& 02126 nfullcols_total=nfullcols_total,& 02127 local_rows=local_rows,& 02128 local_cols=local_cols,& 02129 row_blk_size=row_blk_size,& 02130 col_blk_size=col_blk_size) 02131 02132 !WRITE(*,*)routineN//" sizes",nblkrows_total,& 02133 ! nblkcols_total,& 02134 ! nblkrows_local,& 02135 ! nblkcols_local,& 02136 ! nfullrows_local,& 02137 ! nfullcols_local,& 02138 ! nfullrows_total,& 02139 ! nfullcols_total 02140 02141 02142 rbs => array_data (row_blk_size) 02143 cbs => array_data (col_blk_size) 02144 ALLOCATE (local_row_sizes (nblkrows_total)) 02145 ALLOCATE (local_col_sizes (nblkcols_total)) 02146 local_row_sizes(:) = 0 02147 IF (nblkrows_local .GE. 1) THEN 02148 FORALL (row = 1 : nblkrows_local) 02149 local_row_sizes(local_rows(row)) = rbs(local_rows(row)) 02150 END FORALL 02151 ENDIF 02152 local_col_sizes(:) = 0 02153 IF (nblkcols_local .GE. 1) THEN 02154 FORALL (col = 1 : nblkcols_local) 02155 local_col_sizes(local_cols(col)) = cbs(local_cols(col)) 02156 END FORALL 02157 ENDIF 02158 02159 ALLOCATE (first_row(nblkrows_total),last_row(nblkrows_total)) 02160 ALLOCATE (first_col(nblkcols_total),last_col(nblkcols_total)) 02161 CALL convert_sizes_to_offsets (local_row_sizes, first_row, last_row) 02162 CALL convert_sizes_to_offsets (local_col_sizes, first_col, last_col) 02163 02164 ! Copy the FM data to the block-cyclic DBCSR matrix. This step 02165 ! could be skipped with appropriate DBCSR index manipulation. 02166 fm_block => fm%local_data 02167 fm_block_sp => fm%local_data_sp 02168 02169 CALL cp_dbcsr_work_create (bc_mat, nblks_guess=nblkrows_local*nblkcols_local,& 02170 sizedata_guess=nfullrows_local*nfullcols_local, work_mutable=.FALSE.,& 02171 n=1, error=error) 02172 blk_p = 1 02173 bc_rows: DO row_l = 1, nblkrows_local 02174 row = local_rows (row_l) 02175 row_size = rbs(row) 02176 bc_cols: DO col_l = 1, nblkcols_local 02177 col = local_cols (col_l) 02178 col_size = cbs(col) 02179 nze = row_size*col_size 02180 !WRITE(*,*)routineN//" Adding block",row,col,"size",nze 02181 CALL add_work_coordinate(bc_mat%matrix%m%wms(1), row, col, blk_p, error=dbcsr_error) 02182 IF (fm%use_sp) THEN 02183 !blk_1d_sp => bc_mat%m%wms(1)%data_area%d%r_sp(blk_p:blk_p+nze-1) 02184 blk_1d_sp => dbcsr_get_data_p (bc_mat%matrix%m%wms(1)%data_area,& 02185 coersion=REAL(0.0, KIND=sp), lb=blk_p, ub=blk_p+nze-1) 02186 ELSE 02187 !blk_1d_dp => bc_mat%m%wms(1)%data_area%d%r_dp(blk_p:blk_p+nze-1) 02188 blk_1d_dp => dbcsr_get_data_p (bc_mat%matrix%m%wms(1)%data_area,& 02189 coersion=REAL(0.0, KIND=dp), lb=blk_p, ub=blk_p+nze-1) 02190 ENDIF 02191 CALL cp_assert (nze .EQ. (last_row(row)-first_row(row)+1)*(last_col(col)-first_col(col)+1),& 02192 cp_fatal_level, cp_internal_error, routineN,& 02193 "Block size does not match block row/col sizes") 02194 IF (fm%use_sp) THEN 02195 blk_1d_sp(1:nze) = RESHAPE(& 02196 fm_block_sp(& 02197 first_row(row):last_row(row),first_col(col):last_col(col)& 02198 ), (/ nze /)) 02199 ELSE 02200 blk_1d_dp(1:nze) = RESHAPE(& 02201 fm_block(& 02202 first_row(row):last_row(row),first_col(col):last_col(col)& 02203 ), (/ nze /)) 02204 ENDIF 02205 blk_p = blk_p + nze 02206 ENDDO bc_cols 02207 ENDDO bc_rows 02208 bc_mat%matrix%m%wms(1)%datasize = blk_p - 1 02209 CALL cp_dbcsr_finalize (bc_mat, reshuffle=.FALSE., error=error) 02210 02211 ! Now convert to the desired matrix distribution 02212 IF (PRESENT (alpha)) THEN 02213 CALL stop_program(routineN,moduleN,__LINE__,'no more alpha... clean me') 02214 ELSE 02215 CALL cp_dbcsr_complete_redistribute (bc_mat, matrix,& 02216 keep_sparsity=keep_sparsity, error=error) 02217 ENDIF 02218 CALL cp_dbcsr_release (bc_mat, error=error) 02219 02220 CALL timestop(handle) 02221 CALL dbcsr_error_stop(error_handler, dbcsr_error) 02222 END SUBROUTINE copy_fm_to_dbcsr 02223 02224 ! ***************************************************************************** 02241 SUBROUTINE copy_cfm_to_dbcsr(fm,matrix,keep_sparsity,error) 02242 TYPE(cp_cfm_type), POINTER :: fm 02243 TYPE(cp_dbcsr_type), INTENT(INOUT) :: matrix 02244 LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity 02245 TYPE(cp_error_type), INTENT(inout) :: error 02246 02247 CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_cfm_to_dbcsr', 02248 routineP = moduleN//':'//routineN 02249 02250 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: blk_1d 02251 COMPLEX(KIND=dp), DIMENSION(:, :), 02252 POINTER :: fm_block 02253 INTEGER :: blk_p, col, col_l, col_size, error_handler, group, handle, 02254 nblkcols_local, nblkcols_total, nblkrows_local, nblkrows_total, 02255 ncol_block, ncol_global, nfullcols_local, nfullcols_total, 02256 nfullrows_local, nfullrows_total, nrow_block, nrow_global, nze, row, 02257 row_l, row_size 02258 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, 02259 last_col, last_row, 02260 local_col_sizes, 02261 local_row_sizes 02262 INTEGER, DIMENSION(:), POINTER :: cbs, local_cols, local_rows, 02263 rbs 02264 TYPE(array_i1d_obj) :: col_blk_size, row_blk_size 02265 TYPE(cp_blacs_env_type), POINTER :: context 02266 TYPE(cp_dbcsr_type) :: bc_mat 02267 TYPE(dbcsr_distribution_obj) :: bc_dist 02268 TYPE(dbcsr_error_type) :: dbcsr_error 02269 02270 ! --------------------------------------------------------------------------- 02271 02272 CALL dbcsr_error_set(routineN, error_handler, dbcsr_error) 02273 CALL dbcsr_access_flush (matrix%matrix, error=dbcsr_error) 02274 CALL timeset(routineN,handle) 02275 02276 group = fm%matrix_struct%para_env%group 02277 context => fm%matrix_struct%context 02278 nrow_block = fm%matrix_struct%nrow_block 02279 ncol_block = fm%matrix_struct%ncol_block 02280 nrow_global = fm%matrix_struct%nrow_global 02281 ncol_global = fm%matrix_struct%ncol_global 02282 02283 CALL cp_dbcsr_get_info(matrix,& 02284 nfullrows_total=nfullrows_total,& 02285 nfullcols_total=nfullcols_total) 02286 02287 CALL cp_assert (nrow_global.eq.nfullrows_total,& 02288 cp_fatal_level, cp_caller_error, routineN,& 02289 "FM and DBCSR matrix sizes do not match in rows") 02290 CALL cp_assert (ncol_global.eq.nfullcols_total,& 02291 cp_fatal_level, cp_caller_error, routineN,& 02292 "FM and DBCSR matrix sizes do not match in columns") 02293 02294 ! Create a block-cyclic distribution compatible with the FM matrix. 02295 CALL dbcsr_distribution_init (bc_dist) 02296 CALL dbcsr_create_dist_block_cyclic (bc_dist,& 02297 nfullrows_total, nfullcols_total,& ! Actual full matrix size 02298 nrow_block, ncol_block,& ! BLACS parameters 02299 dbcsr_distribution_mp (cp_dbcsr_distribution (matrix)),& 02300 row_blk_size, col_blk_size) ! block-cyclic row/col sizes 02301 02302 ! Create the block-cyclic DBCSR matrix 02303 CALL cp_dbcsr_init (bc_mat, error) 02304 CALL cp_dbcsr_create (bc_mat, "Block-cyclic "//matrix%matrix%m%name, bc_dist,& 02305 cp_dbcsr_get_matrix_type(matrix), row_blk_size, col_blk_size, 0, 0,& 02306 dbcsr_type_complex_8,error=error) ! type hard coded ! 02307 !call dbcsr_finalize (bc_mat) 02308 CALL dbcsr_distribution_release (bc_dist) 02309 CALL array_release (row_blk_size) 02310 CALL array_release (col_blk_size) 02311 02312 CALL cp_dbcsr_get_info(bc_mat,& 02313 nblkrows_total=nblkrows_total,& 02314 nblkcols_total=nblkcols_total,& 02315 nblkrows_local=nblkrows_local,& 02316 nblkcols_local=nblkcols_local,& 02317 nfullrows_local=nfullrows_local,& 02318 nfullcols_local=nfullcols_local,& 02319 nfullrows_total=nfullrows_total,& 02320 nfullcols_total=nfullcols_total,& 02321 local_rows=local_rows,& 02322 local_cols=local_cols,& 02323 row_blk_size=row_blk_size,& 02324 col_blk_size=col_blk_size) 02325 02326 rbs => array_data (row_blk_size) 02327 cbs => array_data (col_blk_size) 02328 ALLOCATE (local_row_sizes (nblkrows_total)) 02329 ALLOCATE (local_col_sizes (nblkcols_total)) 02330 local_row_sizes(:) = 0 02331 IF (nblkrows_local .GE. 1) THEN 02332 FORALL (row = 1 : nblkrows_local) 02333 local_row_sizes(local_rows(row)) = rbs(local_rows(row)) 02334 END FORALL 02335 ENDIF 02336 local_col_sizes(:) = 0 02337 IF (nblkcols_local .GE. 1) THEN 02338 FORALL (col = 1 : nblkcols_local) 02339 local_col_sizes(local_cols(col)) = cbs(local_cols(col)) 02340 END FORALL 02341 ENDIF 02342 02343 ALLOCATE (first_row(nblkrows_total),last_row(nblkrows_total)) 02344 ALLOCATE (first_col(nblkcols_total),last_col(nblkcols_total)) 02345 CALL convert_sizes_to_offsets (local_row_sizes, first_row, last_row) 02346 CALL convert_sizes_to_offsets (local_col_sizes, first_col, last_col) 02347 02348 ! Copy the FM data to the block-cyclic DBCSR matrix. This step 02349 ! could be skipped with appropriate DBCSR index manipulation. 02350 fm_block => fm%local_data 02351 CALL cp_dbcsr_work_create (bc_mat, nblks_guess=nblkrows_local*nblkcols_local,& 02352 sizedata_guess=nfullrows_local*nfullcols_local, work_mutable=.FALSE.,& 02353 n=1, error=error) 02354 blk_p = 1 02355 bc_rows: DO row_l = 1, nblkrows_local 02356 row = local_rows (row_l) 02357 row_size = rbs(row) 02358 bc_cols: DO col_l = 1, nblkcols_local 02359 col = local_cols (col_l) 02360 col_size = cbs(col) 02361 nze = row_size*col_size 02362 !WRITE(*,*)routineN//" Adding block",row,col,"size",nze 02363 CALL add_work_coordinate(bc_mat%matrix%m%wms(1), row, col, blk_p, error=dbcsr_error) 02364 !blk_1d => bc_mat%m%wms(1)%data_area%d%c_dp(blk_p:blk_p+nze-1) 02365 blk_1d => dbcsr_get_data_p (bc_mat%matrix%m%wms(1)%data_area,& 02366 coersion=CMPLX(0.0, KIND=dp), lb=blk_p, ub=blk_p+nze-1) 02367 CALL cp_assert (nze .EQ. (last_row(row)-first_row(row)+1)*(last_col(col)-first_col(col)+1),& 02368 cp_fatal_level, cp_internal_error, routineN,& 02369 "Block size does not match block row/col sizes") 02370 blk_1d(1:nze) = RESHAPE(& 02371 fm_block(& 02372 first_row(row):last_row(row),first_col(col):last_col(col)& 02373 ), (/ nze /)) 02374 blk_p = blk_p + nze 02375 ENDDO bc_cols 02376 ENDDO bc_rows 02377 bc_mat%matrix%m%wms(1)%datasize = blk_p - 1 02378 CALL cp_dbcsr_finalize (bc_mat, reshuffle=.FALSE., error=error) 02379 02380 ! Now convert to the desired matrix distribution 02381 CALL cp_dbcsr_complete_redistribute (bc_mat, matrix, keep_sparsity=keep_sparsity,& 02382 error=error) 02383 CALL cp_dbcsr_release (bc_mat, error=error) 02384 02385 CALL timestop(handle) 02386 CALL dbcsr_error_stop(error_handler, dbcsr_error) 02387 END SUBROUTINE copy_cfm_to_dbcsr 02388 02389 ! ***************************************************************************** 02394 SUBROUTINE copy_dbcsr_to_fm(matrix, fm, error) 02395 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix 02396 TYPE(cp_fm_type), POINTER :: fm 02397 TYPE(cp_error_type), INTENT(INOUT) :: error 02398 02399 CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_dbcsr_to_fm', 02400 routineP = moduleN//':'//routineN 02401 LOGICAL, PARAMETER :: dbg = .FALSE. 02402 02403 INTEGER :: col, error_handle, group, handle, mypcol, mype, myprow, 02404 nblkcols_local, nblkcols_total, nblkrows_local, nblkrows_total, 02405 ncol_block, ncol_global, nfullcols_total, nfullrows_total, npcol, npe, 02406 nprow, nrow_block, nrow_global, row 02407 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, 02408 last_col, last_row, 02409 local_col_sizes, 02410 local_row_sizes 02411 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, local_cols, 02412 local_rows, row_blk_sizes 02413 REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk_2d, fm_block 02414 REAL(KIND=sp), DIMENSION(:, :), POINTER :: fm_block_sp 02415 TYPE(array_i1d_obj) :: cbs, rbs 02416 TYPE(cp_blacs_env_type), POINTER :: context 02417 TYPE(cp_dbcsr_iterator) :: iter 02418 TYPE(cp_dbcsr_type) :: bc_mat 02419 TYPE(dbcsr_distribution_obj) :: bc_dist 02420 TYPE(dbcsr_error_type) :: dbcsr_error 02421 02422 ! --------------------------------------------------------------------------- 02423 02424 CALL timeset(routineN,handle) 02425 CALL dbcsr_error_set (routineN, error_handle, error=dbcsr_error) 02426 CALL dbcsr_access_flush (matrix%matrix, error=dbcsr_error) 02427 02428 ! info about the full matrix 02429 group = fm%matrix_struct%para_env%group 02430 context => fm%matrix_struct%context 02431 mype=context%my_pid 02432 npe=context%n_pid 02433 myprow=context%mepos(1) 02434 mypcol=context%mepos(2) 02435 nprow=context%num_pe(1) 02436 npcol=context%num_pe(2) 02437 nrow_block = fm%matrix_struct%nrow_block 02438 ncol_block = fm%matrix_struct%ncol_block 02439 nrow_global = fm%matrix_struct%nrow_global 02440 ncol_global = fm%matrix_struct%ncol_global 02441 02442 CALL cp_dbcsr_get_info(matrix,& 02443 nfullrows_total=nfullrows_total,& 02444 nfullcols_total=nfullcols_total) 02445 02446 ! Convert DBCSR to a block-cyclic one 02447 CALL dbcsr_distribution_init (bc_dist) 02448 CALL dbcsr_create_dist_block_cyclic (bc_dist,& 02449 nfullrows_total, nfullcols_total,& 02450 nrow_block, ncol_block,& 02451 dbcsr_distribution_mp (cp_dbcsr_distribution(matrix)),& 02452 rbs, cbs) 02453 02454 CALL cp_dbcsr_init (bc_mat, error) 02455 CALL cp_dbcsr_create (bc_mat, "Block-cyclic"//matrix%matrix%m%name, bc_dist,& 02456 dbcsr_type_no_symmetry, rbs, cbs, 0, 0, error=error) 02457 CALL dbcsr_distribution_release (bc_dist) 02458 CALL array_release (rbs) 02459 CALL array_release (cbs) 02460 CALL cp_dbcsr_complete_redistribute (matrix, bc_mat, error=error) 02461 02462 ! Find the local extents of the local blocked rows so that index lookups 02463 ! into the FM matrix work correctly. 02464 row_blk_sizes => array_data (rbs) 02465 col_blk_sizes => array_data (cbs) 02466 local_rows => array_data (dbcsr_distribution_local_rows (bc_dist)) 02467 local_cols => array_data (dbcsr_distribution_local_cols (bc_dist)) 02468 ALLOCATE (local_row_sizes (dbcsr_distribution_nrows (bc_dist))) 02469 ALLOCATE (local_col_sizes (dbcsr_distribution_ncols (bc_dist))) 02470 nblkrows_local = dbcsr_distribution_nlocal_rows (bc_dist) 02471 nblkcols_local = dbcsr_distribution_nlocal_cols (bc_dist) 02472 local_row_sizes(:) = 0 02473 IF (nblkrows_local .GE. 1) THEN 02474 FORALL (row = 1 : nblkrows_local) 02475 local_row_sizes(local_rows(row)) = row_blk_sizes(local_rows(row)) 02476 END FORALL 02477 ENDIF 02478 local_col_sizes(:) = 0 02479 IF (nblkcols_local .GE. 1) THEN 02480 FORALL (col = 1 : nblkcols_local) 02481 local_col_sizes(local_cols(col)) = col_blk_sizes(local_cols(col)) 02482 END FORALL 02483 ENDIF 02484 nblkrows_total = dbcsr_distribution_nrows (bc_dist) 02485 nblkcols_total = dbcsr_distribution_ncols (bc_dist) 02486 ALLOCATE (first_row(nblkrows_total),last_row(nblkrows_total)) 02487 ALLOCATE (first_col(nblkcols_total),last_col(nblkcols_total)) 02488 CALL convert_sizes_to_offsets (local_row_sizes, first_row, last_row) 02489 CALL convert_sizes_to_offsets (local_col_sizes, first_col, last_col) 02490 ! 02491 ! Now copy data to the FM matrix 02492 fm_block => fm%local_data 02493 fm_block_sp => fm%local_data_sp 02494 IF(fm%use_sp) THEN 02495 fm_block_sp=0.0_sp 02496 ELSE 02497 fm_block=0.0_dp 02498 ENDIF 02499 02500 IF(dbg) THEN 02501 WRITE(*,*)routineN//" FM data size is", UBOUND(fm_block) 02502 WRITE(*,*)routineN//" dbcsr data size is", cp_dbcsr_get_data_size(bc_mat) 02503 WRITE(*,*)routineN//" FM block sizes are",nrow_block,'/',nfullrows_total 02504 WRITE(*,*)routineN//" FM block sizes are",ncol_block,'/',nfullcols_total 02505 WRITE(*,*)routineN//" dbcsr row sizes are",bc_mat%matrix%m%row_blk_size%low%data 02506 WRITE(*,*)routineN//" dbcsr col sizes are",bc_mat%matrix%m%col_blk_size%low%data 02507 ENDIF 02508 ! 02509 CALL cp_dbcsr_iterator_start(iter, bc_mat) 02510 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 02511 CALL cp_dbcsr_iterator_next_block (iter, row, col, blk_2d) 02512 ! Convert absolute coordinates to FM-local coordinates 02513 IF(fm%use_sp) THEN 02514 fm_block_sp(first_row(row):last_row(row),first_col(col):last_col(col))& 02515 = REAL(blk_2d(:,:),sp) 02516 ELSE 02517 IF (dbg) THEN 02518 WRITE(*,*)routineN//" blk2d size",UBOUND(blk_2d) 02519 WRITE(*,*)routineN//" want to set coor.",row,col 02520 WRITE(*,*)routineN//" local extents",& 02521 first_row(row),last_row(row),first_col(col),last_col(col) 02522 ENDIF 02523 fm_block(first_row(row):last_row(row),first_col(col):last_col(col))& 02524 = blk_2d(:,:) 02525 ENDIF 02526 ENDDO 02527 CALL cp_dbcsr_iterator_stop(iter) 02528 02529 CALL cp_dbcsr_release (bc_mat, error=error) 02530 CALL dbcsr_error_stop (error_handle, error=dbcsr_error) 02531 CALL timestop(handle) 02532 END SUBROUTINE copy_dbcsr_to_fm 02533 02534 ! ***************************************************************************** 02539 SUBROUTINE copy_dbcsr_to_cfm(matrix, fm, error) 02540 TYPE(cp_dbcsr_type), INTENT(IN) :: matrix 02541 TYPE(cp_cfm_type), POINTER :: fm 02542 TYPE(cp_error_type), INTENT(INOUT) :: error 02543 02544 CHARACTER(LEN=*), PARAMETER :: routineN = 'copy_dbcsr_to_cfm', 02545 routineP = moduleN//':'//routineN 02546 LOGICAL, PARAMETER :: dbg = .FALSE. 02547 02548 COMPLEX(KIND=dp), DIMENSION(:, :), 02549 POINTER :: blk_2d, fm_block 02550 INTEGER :: col, error_handle, group, handle, mypcol, mype, myprow, 02551 nblkcols_local, nblkcols_total, nblkrows_local, nblkrows_total, 02552 ncol_block, ncol_global, nfullcols_total, nfullrows_total, npcol, npe, 02553 nprow, nrow_block, nrow_global, row 02554 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_col, first_row, 02555 last_col, last_row, 02556 local_col_sizes, 02557 local_row_sizes 02558 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, local_cols, 02559 local_rows, row_blk_sizes 02560 TYPE(array_i1d_obj) :: cbs, rbs 02561 TYPE(cp_blacs_env_type), POINTER :: context 02562 TYPE(cp_dbcsr_iterator) :: iter 02563 TYPE(cp_dbcsr_type) :: bc_mat 02564 TYPE(dbcsr_distribution_obj) :: bc_dist 02565 TYPE(dbcsr_error_type) :: dbcsr_error 02566 02567 ! --------------------------------------------------------------------------- 02568 02569 CALL timeset(routineN,handle) 02570 CALL dbcsr_error_set (routineN, error_handle, error=dbcsr_error) 02571 CALL dbcsr_access_flush (matrix%matrix, error=dbcsr_error) 02572 02573 ! info about the full matrix 02574 group = fm%matrix_struct%para_env%group 02575 context => fm%matrix_struct%context 02576 mype=context%my_pid 02577 npe=context%n_pid 02578 myprow=context%mepos(1) 02579 mypcol=context%mepos(2) 02580 nprow=context%num_pe(1) 02581 npcol=context%num_pe(2) 02582 nrow_block = fm%matrix_struct%nrow_block 02583 ncol_block = fm%matrix_struct%ncol_block 02584 nrow_global = fm%matrix_struct%nrow_global 02585 ncol_global = fm%matrix_struct%ncol_global 02586 02587 CALL cp_dbcsr_get_info(matrix,& 02588 nfullrows_total=nfullrows_total,& 02589 nfullcols_total=nfullcols_total) 02590 02591 ! Convert DBCSR to a block-cyclic one 02592 CALL dbcsr_distribution_init (bc_dist) 02593 CALL dbcsr_create_dist_block_cyclic (bc_dist,& 02594 nfullrows_total, nfullcols_total,& 02595 nrow_block, ncol_block,& 02596 dbcsr_distribution_mp (cp_dbcsr_distribution(matrix)),& 02597 rbs, cbs) 02598 02599 CALL cp_dbcsr_init (bc_mat, error) 02600 CALL cp_dbcsr_create (bc_mat, "Block-cyclic"//matrix%matrix%m%name, bc_dist,& 02601 dbcsr_type_no_symmetry, rbs, cbs, 0, 0, cp_dbcsr_get_data_type(matrix),& 02602 error=error) 02603 CALL dbcsr_distribution_release (bc_dist) 02604 CALL array_release (rbs) 02605 CALL array_release (cbs) 02606 CALL cp_dbcsr_complete_redistribute (matrix, bc_mat, error=error) 02607 02608 ! Find the local extents of the local blocked rows so that index lookups 02609 ! into the FM matrix work correctly. 02610 row_blk_sizes => array_data (rbs) 02611 col_blk_sizes => array_data (cbs) 02612 local_rows => array_data (dbcsr_distribution_local_rows (bc_dist)) 02613 local_cols => array_data (dbcsr_distribution_local_cols (bc_dist)) 02614 ALLOCATE (local_row_sizes (dbcsr_distribution_nrows (bc_dist))) 02615 ALLOCATE (local_col_sizes (dbcsr_distribution_ncols (bc_dist))) 02616 nblkrows_local = dbcsr_distribution_nlocal_rows (bc_dist) 02617 nblkcols_local = dbcsr_distribution_nlocal_cols (bc_dist) 02618 local_row_sizes(:) = 0 02619 IF (nblkrows_local .GE. 1) THEN 02620 FORALL (row = 1 : nblkrows_local) 02621 local_row_sizes(local_rows(row)) = row_blk_sizes(local_rows(row)) 02622 END FORALL 02623 ENDIF 02624 local_col_sizes(:) = 0 02625 IF (nblkcols_local .GE. 1) THEN 02626 FORALL (col = 1 : nblkcols_local) 02627 local_col_sizes(local_cols(col)) = col_blk_sizes(local_cols(col)) 02628 END FORALL 02629 ENDIF 02630 nblkrows_total = dbcsr_distribution_nrows (bc_dist) 02631 nblkcols_total = dbcsr_distribution_ncols (bc_dist) 02632 ALLOCATE (first_row(nblkrows_total),last_row(nblkrows_total)) 02633 ALLOCATE (first_col(nblkcols_total),last_col(nblkcols_total)) 02634 CALL convert_sizes_to_offsets (local_row_sizes, first_row, last_row) 02635 CALL convert_sizes_to_offsets (local_col_sizes, first_col, last_col) 02636 ! 02637 ! Now copy data to the FM matrix 02638 fm_block => fm%local_data 02639 fm_block=(0.0_dp,0.0_dp) 02640 02641 IF(dbg) THEN 02642 WRITE(*,*)routineN//" FM data size is", UBOUND(fm_block) 02643 WRITE(*,*)routineN//" dbcsr data size is", cp_dbcsr_get_data_size(bc_mat) 02644 WRITE(*,*)routineN//" FM block sizes are",nrow_block,'/',nfullrows_total 02645 WRITE(*,*)routineN//" FM block sizes are",ncol_block,'/',nfullcols_total 02646 WRITE(*,*)routineN//" dbcsr row sizes are",bc_mat%matrix%m%row_blk_size%low%data 02647 WRITE(*,*)routineN//" dbcsr col sizes are",bc_mat%matrix%m%col_blk_size%low%data 02648 ENDIF 02649 ! 02650 CALL cp_dbcsr_iterator_start(iter, bc_mat) 02651 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 02652 CALL cp_dbcsr_iterator_next_block (iter, row, col, blk_2d) 02653 ! Convert absolute coordinates to FM-local coordinates 02654 IF (dbg) THEN 02655 WRITE(*,*)routineN//" blk2d size",UBOUND(blk_2d) 02656 WRITE(*,*)routineN//" want to set coor.",row,col 02657 WRITE(*,*)routineN//" local extents",& 02658 first_row(row),last_row(row),first_col(col),last_col(col) 02659 ENDIF 02660 fm_block(first_row(row):last_row(row),first_col(col):last_col(col))& 02661 = blk_2d(:,:) 02662 ENDDO 02663 CALL cp_dbcsr_iterator_stop(iter) 02664 02665 CALL cp_dbcsr_release (bc_mat, error=error) 02666 CALL dbcsr_error_stop (error_handle, error=dbcsr_error) 02667 CALL timestop(handle) 02668 END SUBROUTINE copy_dbcsr_to_cfm 02669 02670 END MODULE cp_dbcsr_operations
1.7.3