|
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 00016 MODULE dbcsr_operations 00017 USE array_types, ONLY: array_data,& 00018 array_equality,& 00019 array_get,& 00020 array_hold,& 00021 array_i1d_obj,& 00022 array_release 00023 USE dbcsr_blas_operations, ONLY: dbcsr_blas_copy,& 00024 dbcsr_blas_gemm 00025 USE dbcsr_block_access, ONLY: dbcsr_access_flush,& 00026 dbcsr_get_block_p,& 00027 dbcsr_put_block,& 00028 dbcsr_remove_block,& 00029 dbcsr_reserve_blocks 00030 USE dbcsr_block_operations, ONLY: block_add_on_diag,& 00031 dbcsr_block_conjg,& 00032 dbcsr_block_partial_copy,& 00033 dbcsr_block_real_neg,& 00034 dbcsr_block_scale,& 00035 dbcsr_data_clear,& 00036 get_block2d_diagonal,& 00037 set_block2d_diagonal 00038 USE dbcsr_config, ONLY: dbcsr_init_conf,& 00039 has_cuda,& 00040 is_configured,& 00041 mm_async,& 00042 mm_driver,& 00043 mm_driver_plasma 00044 USE dbcsr_cuda_device, ONLY: dbcsr_cuda_get_n_devices,& 00045 dbcsr_cuda_init 00046 USE dbcsr_data_methods, ONLY: & 00047 dbcsr_data_clear_pointer, dbcsr_data_ensure_size, dbcsr_data_get_size, & 00048 dbcsr_data_get_type, dbcsr_data_init, dbcsr_data_new, & 00049 dbcsr_data_release, dbcsr_data_set_pointer, & 00050 dbcsr_data_set_size_referenced, dbcsr_get_data, dbcsr_scalar, & 00051 dbcsr_scalar_are_equal, dbcsr_scalar_fill_all, dbcsr_scalar_get_type, & 00052 dbcsr_scalar_one, dbcsr_scalar_set_type, dbcsr_scalar_zero, & 00053 dbcsr_type_1d_to_2d 00054 USE dbcsr_data_operations, ONLY: dbcsr_data_convert,& 00055 dbcsr_data_copyall 00056 USE dbcsr_dist_operations, ONLY: & 00057 checker_square_proc, checker_tr, create_bl_distribution, & 00058 dbcsr_create_image_dist, dbcsr_find_column, & 00059 dbcsr_get_stored_coordinates, dbcsr_make_dists_dense, & 00060 dbcsr_reset_locals, dbcsr_reset_vlocals, make_sizes_dense 00061 USE dbcsr_error_handling 00062 USE dbcsr_index_operations, ONLY: dbcsr_expand_row_index,& 00063 dbcsr_index_checksum,& 00064 dbcsr_index_compact,& 00065 dbcsr_make_index_canonical,& 00066 dbcsr_make_index_list,& 00067 dbcsr_make_index_local_row,& 00068 dbcsr_repoint_index 00069 USE dbcsr_io, ONLY: dbcsr_print 00070 USE dbcsr_iterator_operations, ONLY: dbcsr_iterator_blocks_left,& 00071 dbcsr_iterator_next_block,& 00072 dbcsr_iterator_start,& 00073 dbcsr_iterator_stop 00074 USE dbcsr_kinds, ONLY: & 00075 dp, int_1_size, int_2_size, int_4_size, int_8, int_8_size, real_4, & 00076 real_8, sp 00077 USE dbcsr_machine, ONLY: default_output_unit 00078 USE dbcsr_message_passing, ONLY: dmp_max,& 00079 mp_allgather,& 00080 mp_bcast,& 00081 mp_environ,& 00082 mp_recv,& 00083 mp_send,& 00084 mp_sum 00085 USE dbcsr_methods, ONLY: & 00086 dbcsr_col_block_offsets, dbcsr_col_block_sizes, dbcsr_destroy_array, & 00087 dbcsr_distribution, dbcsr_distribution_col_dist, & 00088 dbcsr_distribution_has_threads, dbcsr_distribution_hold, & 00089 dbcsr_distribution_make_threads, dbcsr_distribution_mp, & 00090 dbcsr_distribution_ncols, dbcsr_distribution_new, & 00091 dbcsr_distribution_no_threads, dbcsr_distribution_nrows, & 00092 dbcsr_distribution_release, dbcsr_distribution_row_dist, & 00093 dbcsr_get_data_size, dbcsr_get_data_size_referenced, & 00094 dbcsr_get_data_type, dbcsr_get_index_memory_type, dbcsr_get_info, & 00095 dbcsr_get_matrix_type, dbcsr_get_num_blocks, & 00096 dbcsr_get_replication_type, dbcsr_has_symmetry, dbcsr_image_dist_hold, & 00097 dbcsr_image_dist_release, dbcsr_init, dbcsr_is_initialized, & 00098 dbcsr_max_col_size, dbcsr_max_row_size, dbcsr_may_be_dense, & 00099 dbcsr_mp_grid_setup, dbcsr_mp_group, dbcsr_mp_has_subgroups, & 00100 dbcsr_mp_my_col_group, dbcsr_mp_my_row_group, dbcsr_mp_mynode, & 00101 dbcsr_mp_mypcol, dbcsr_mp_myprow, dbcsr_mp_npcols, dbcsr_mp_nprows, & 00102 dbcsr_mp_numnodes, dbcsr_mp_pgrid, dbcsr_name, dbcsr_nblkcols_total, & 00103 dbcsr_nblkrows_total, dbcsr_nfullcols_total, dbcsr_nfullrows_total, & 00104 dbcsr_release, dbcsr_release_locals, dbcsr_row_block_offsets, & 00105 dbcsr_row_block_sizes, dbcsr_switch_data_area, dbcsr_valid_index 00106 USE dbcsr_mm_cannon, ONLY: dbcsr_mm_cannon_lib_finalize,& 00107 dbcsr_mm_cannon_lib_init,& 00108 dbcsr_mm_cannon_multiply 00109 USE dbcsr_mp_operations, ONLY: dbcsr_recv_any,& 00110 dbcsr_send_any 00111 USE dbcsr_plasma_interface, ONLY: dbcsr_plasma_finalize,& 00112 dbcsr_plasma_init 00113 USE dbcsr_ptr_util, ONLY: ensure_array_size,& 00114 pointer_view 00115 USE dbcsr_toollib, ONLY: ceil_log2,& 00116 swap,& 00117 uppercase 00118 USE dbcsr_transformations, ONLY: dbcsr_make_dense,& 00119 dbcsr_make_images,& 00120 dbcsr_make_images_dense,& 00121 dbcsr_make_undense,& 00122 dbcsr_make_untransposed_blocks,& 00123 dbcsr_new_transposed 00124 USE dbcsr_types, ONLY: & 00125 dbcsr_2d_array_type, dbcsr_conjugate_transpose, dbcsr_data_obj, & 00126 dbcsr_distribution_obj, dbcsr_filter_frobenius, & 00127 dbcsr_imagedistribution_obj, dbcsr_iterator, dbcsr_memory_MPI, & 00128 dbcsr_memory_default, dbcsr_mp_obj, dbcsr_no_transpose, & 00129 dbcsr_norm_column, dbcsr_norm_frobenius, dbcsr_norm_gershgorin, & 00130 dbcsr_norm_maxabsnorm, dbcsr_obj, dbcsr_repl_col, dbcsr_repl_full, & 00131 dbcsr_repl_none, dbcsr_repl_row, dbcsr_scalar_type, dbcsr_transpose, & 00132 dbcsr_type, dbcsr_type_antihermitian, dbcsr_type_antisymmetric, & 00133 dbcsr_type_complex_4, dbcsr_type_complex_8, dbcsr_type_hermitian, & 00134 dbcsr_type_no_symmetry, dbcsr_type_real_4, dbcsr_type_real_8, & 00135 dbcsr_type_symmetric 00136 USE dbcsr_util, ONLY: dbcsr_checksum,& 00137 dbcsr_verify_matrix,& 00138 find_block_of_element 00139 USE dbcsr_work_operations, ONLY: dbcsr_add_wm_from_matrix,& 00140 dbcsr_create,& 00141 dbcsr_finalize,& 00142 dbcsr_work_create 00143 USE ma_affinity, ONLY: ma_set_gpu_affinity 00144 USE machine_architecture_types, ONLY: has_ma 00145 00146 !$ USE OMP_LIB 00147 00148 IMPLICIT NONE 00149 00150 PRIVATE 00151 00152 00153 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbcsr_operations' 00154 00155 REAL, PARAMETER :: default_resize_factor = 1.618034 00156 00157 ! prettify protection 00158 CHARACTER, PARAMETER :: xa=dbcsr_type_hermitian, xb=dbcsr_type_antihermitian 00159 00160 PUBLIC :: dbcsr_init_lib, dbcsr_finalize_lib 00161 PUBLIC :: dbcsr_multiply,& 00162 dbcsr_trace, dbcsr_add_on_diag,& 00163 dbcsr_set, dbcsr_scale, dbcsr_scale_mat, dbcsr_add, dbcsr_copy,& 00164 dbcsr_copy_submatrix, dbcsr_copy_into_existing,& 00165 dbcsr_get_diag, dbcsr_set_diag, & 00166 dbcsr_get_block_diag, dbcsr_hadamard_product, & 00167 dbcsr_filter, dbcsr_scale_by_vector,& 00168 dbcsr_replace_blocks, dbcsr_conjg,& 00169 dbcsr_btriu, dbcsr_triu, dbcsr_tril,& 00170 dbcsr_symmetrize_block_diag, dbcsr_copy_columns,& 00171 dbcsr_init_random, dbcsr_lanczos, dbcsr_block_in_limits 00172 PUBLIC :: dbcsr_sum_replicated 00173 PUBLIC :: dbcsr_norm, & 00174 dbcsr_gershgorin_norm, dbcsr_maxabs, dbcsr_frobenius_norm 00175 00176 ! The interfaces for the generic routines found in the generated 00177 ! generic files. 00178 00179 INTERFACE dbcsr_multiply 00180 MODULE PROCEDURE dbcsr_multiply_anytype 00181 MODULE PROCEDURE dbcsr_multiply_s, dbcsr_multiply_d,& 00182 dbcsr_multiply_c, dbcsr_multiply_z 00183 END INTERFACE 00184 00185 INTERFACE dbcsr_trace 00186 MODULE PROCEDURE dbcsr_trace_a_any 00187 MODULE PROCEDURE dbcsr_trace_a_s, dbcsr_trace_a_d,& 00188 dbcsr_trace_a_c, dbcsr_trace_a_z 00189 MODULE PROCEDURE dbcsr_trace_ab_s, dbcsr_trace_a_b_d,& 00190 dbcsr_trace_ab_c, dbcsr_trace_ab_z 00191 END INTERFACE 00192 00193 INTERFACE dbcsr_scale 00194 MODULE PROCEDURE dbcsr_scale_anytype 00195 MODULE PROCEDURE dbcsr_scale_s, dbcsr_scale_d,& 00196 dbcsr_scale_c, dbcsr_scale_z 00197 END INTERFACE 00198 00199 INTERFACE dbcsr_scale_mat 00200 MODULE PROCEDURE dbcsr_scale_mat_any 00201 MODULE PROCEDURE dbcsr_scale_s_m, dbcsr_scale_d_m,& 00202 dbcsr_scale_c_m, dbcsr_scale_z_m 00203 END INTERFACE 00204 00205 INTERFACE dbcsr_scale_by_vector 00206 MODULE PROCEDURE dbcsr_scale_by_vector_anytype 00207 MODULE PROCEDURE dbcsr_scale_by_vector_s, dbcsr_scale_by_vector_d,& 00208 dbcsr_scale_by_vector_c, dbcsr_scale_by_vector_z 00209 END INTERFACE 00210 00211 INTERFACE dbcsr_set 00212 MODULE PROCEDURE dbcsr_set_anytype 00213 MODULE PROCEDURE dbcsr_set_s, dbcsr_set_d, dbcsr_set_c, dbcsr_set_z 00214 END INTERFACE 00215 00216 INTERFACE dbcsr_add 00217 MODULE PROCEDURE dbcsr_add_anytype 00218 MODULE PROCEDURE dbcsr_add_s, dbcsr_add_d,& 00219 dbcsr_add_c, dbcsr_add_z 00220 END INTERFACE 00221 00222 INTERFACE dbcsr_add_on_diag 00223 MODULE PROCEDURE dbcsr_add_on_diag_anytype 00224 MODULE PROCEDURE dbcsr_add_on_diag_s, dbcsr_add_on_diag_d,& 00225 dbcsr_add_on_diag_c, dbcsr_add_on_diag_z 00226 END INTERFACE 00227 00228 INTERFACE dbcsr_filter 00229 MODULE PROCEDURE dbcsr_filter_anytype 00230 MODULE PROCEDURE dbcsr_filter_s, dbcsr_filter_d,& 00231 dbcsr_filter_c, dbcsr_filter_z 00232 END INTERFACE 00233 00234 INTERFACE dbcsr_get_diag 00235 MODULE PROCEDURE dbcsr_get_diag_anytype 00236 MODULE PROCEDURE dbcsr_get_diag_s, dbcsr_get_diag_d,& 00237 dbcsr_get_diag_c, dbcsr_get_diag_z 00238 END INTERFACE 00239 00240 INTERFACE dbcsr_set_diag 00241 MODULE PROCEDURE dbcsr_set_diag_anytype 00242 MODULE PROCEDURE dbcsr_set_diag_s, dbcsr_set_diag_d,& 00243 dbcsr_set_diag_c, dbcsr_set_diag_z 00244 END INTERFACE 00245 00246 INTERFACE dbcsr_norm 00247 MODULE PROCEDURE dbcsr_norm_anytype 00248 MODULE PROCEDURE dbcsr_norm_r4_scal 00249 MODULE PROCEDURE dbcsr_norm_r4_vec, dbcsr_norm_r8_vec 00250 END INTERFACE 00251 INTERFACE dbcsr_gershgorin_norm 00252 MODULE PROCEDURE dbcsr_gershgorin_norm_r8 00253 END INTERFACE 00254 INTERFACE dbcsr_maxabs 00255 MODULE PROCEDURE dbcsr_maxabs_r8 00256 END INTERFACE 00257 INTERFACE dbcsr_frobenius_norm 00258 MODULE PROCEDURE dbcsr_frobenius_norm_r8 00259 END INTERFACE 00260 00261 INTERFACE dbcsr_lanczos 00262 MODULE PROCEDURE dbcsr_lanczos_extremal_eig 00263 END INTERFACE 00264 00265 LOGICAL, PARAMETER :: debug_mod = .FALSE. 00266 LOGICAL, PARAMETER :: careful_mod = .FALSE. 00267 00268 #define temp_transpose(v, r, c) RESHAPE(TRANSPOSE(RESHAPE(v,(/r,c/))),(/r*c/)) 00269 00270 INTEGER, PARAMETER, PRIVATE :: rpslot_owner = 1 00271 INTEGER, PARAMETER, PRIVATE :: rpslot_addblks = 2 00272 INTEGER, PARAMETER, PRIVATE :: rpslot_addoffset = 3 00273 INTEGER, PARAMETER, PRIVATE :: rpslot_oldblks = 4 00274 INTEGER, PARAMETER, PRIVATE :: rpslot_oldoffset = 5 00275 INTEGER, PARAMETER, PRIVATE :: rpslot_totaloffset = 6 00276 INTEGER, PARAMETER, PRIVATE :: rpnslots = 6 00277 00278 00279 CONTAINS 00280 00281 00282 ! ***************************************************************************** 00339 SUBROUTINE dbcsr_multiply_anytype(transa, transb,& 00340 alpha, matrix_a, matrix_b, beta, matrix_c,& 00341 first_row, last_row, first_column, last_column, first_k, last_k,& 00342 retain_sparsity, filter_eps,& 00343 left_set, right_set, error, flop) 00344 00345 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 00346 TYPE(dbcsr_scalar_type), INTENT(IN) :: alpha 00347 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 00348 TYPE(dbcsr_scalar_type), INTENT(IN) :: beta 00349 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 00350 INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, 00351 first_column, last_column, 00352 first_k, last_k 00353 LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity 00354 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: filter_eps 00355 TYPE(dbcsr_2d_array_type), OPTIONAL, 00356 POINTER :: left_set, right_set 00357 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 00358 INTEGER(KIND=int_8), INTENT(OUT), 00359 OPTIONAL :: flop 00360 00361 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_multiply_anytype', 00362 routineP = moduleN//':'//routineN 00363 LOGICAL, PARAMETER :: dbg = .FALSE., 00364 use_list_indexing = .TRUE., 00365 use_local_indexing = .TRUE. 00366 REAL(real_8), PARAMETER :: make_dense_occ_thresh = 1.0_dp 00367 00368 CHARACTER :: transa_l, transb_l 00369 INTEGER :: comm, error_handler, f_col, 00370 f_k, f_row, i, l_col, l_k, 00371 l_row, numnodes, output_unit 00372 INTEGER(KIND=int_8) :: my_flop 00373 LOGICAL :: ab_dense, keep_product_data, keep_sparsity, left_set_given, 00374 new_left, new_right, plasma_is_set, product_reindex, release_tdist, 00375 right_set_given, use_dense_mult, use_plasma 00376 REAL(KIND=dp) :: cs 00377 TYPE(array_i1d_obj) :: dense_col_sizes, dense_k_sizes, dense_row_sizes, 00378 k_vmap, m_map, n_map, old_product_col_blk_offsets, 00379 old_product_col_blk_sizes, old_product_row_blk_offsets, 00380 old_product_row_blk_sizes 00381 TYPE(dbcsr_2d_array_type), POINTER :: m2s_left, m2s_right 00382 TYPE(dbcsr_distribution_obj) :: dense_product_distribution, 00383 old_product_distribution 00384 TYPE(dbcsr_imagedistribution_obj) :: dense_rdist_left, 00385 dense_rdist_right, 00386 rdist_left, rdist_right 00387 TYPE(dbcsr_obj) :: dense_template_left, dense_template_right, 00388 matrix_left, matrix_right, matrix_tmp, product_matrix 00389 TYPE(dbcsr_scalar_type) :: eps_any 00390 00391 CALL dbcsr_error_set(routineN, error_handler, error) 00392 ! 00393 ! check parameters 00394 transa_l = transa 00395 transb_l = transb 00396 CALL uppercase(transa_l) 00397 CALL uppercase(transb_l) 00398 CALL dbcsr_assert(transa_l.EQ.dbcsr_no_transpose.OR.& 00399 transa_l.EQ.dbcsr_transpose.OR.& 00400 transa_l.EQ.dbcsr_conjugate_transpose,& 00401 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00402 "Invalid transa_l = "//transa_l, __LINE__, error) 00403 CALL dbcsr_assert(transb_l.EQ.dbcsr_no_transpose.OR.& 00404 transb_l.EQ.dbcsr_transpose.OR.& 00405 transb_l.EQ.dbcsr_conjugate_transpose,& 00406 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00407 "Invalid transb_l = "//transb_l, __LINE__, error) 00408 00409 IF (dbg) THEN 00410 WRITE(*,*)'========== MULTIPLICATION ========================' 00411 CALL dbcsr_verify_matrix (matrix_a, error=error) 00412 CALL dbcsr_verify_matrix (matrix_b, error=error) 00413 CALL dbcsr_verify_matrix (matrix_c, error=error) 00414 WRITE(*,*)routineN//" ABC checksums",& 00415 dbcsr_checksum(matrix_a, error=error),& 00416 dbcsr_checksum(matrix_b, error=error),& 00417 dbcsr_checksum(matrix_c, error=error) 00418 IF (dbg) THEN 00419 CALL dbcsr_print (matrix_a, nodata=.TRUE., error=error) 00420 CALL dbcsr_print (matrix_b, nodata=.TRUE., error=error) 00421 CALL dbcsr_print (matrix_c, nodata=.TRUE., error=error) 00422 ENDIF 00423 ENDIF 00424 ! 00425 CALL dbcsr_access_flush (matrix_a, error=error) 00426 CALL dbcsr_access_flush (matrix_b, error=error) 00427 CALL dbcsr_access_flush (matrix_c, error=error) 00428 ! 00429 ! transpose/conjg left and/or right matrices if needed 00430 SELECT CASE(transa_l) 00431 CASE(dbcsr_no_transpose) 00432 matrix_left = matrix_a 00433 new_left = .FALSE. 00434 CASE(dbcsr_transpose) 00435 CALL dbcsr_init(matrix_left) 00436 IF(dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric) THEN 00437 ! 00438 ! For antisymmetric matrix, we need to do a hard copy 00439 ! shallow_data_copy=.TRUE. doesnt handle properly antisymm matrices 00440 CALL dbcsr_new_transposed (matrix_left, matrix_a,& 00441 shallow_data_copy=.FALSE., redistribute=.FALSE., & 00442 transpose_distribution=.FALSE., error=error) 00443 ELSE 00444 CALL dbcsr_new_transposed (matrix_left, matrix_a,& 00445 shallow_data_copy=.TRUE., redistribute=.FALSE.,& 00446 transpose_distribution=.FALSE., error=error) 00447 ENDIF 00448 new_left = .TRUE. 00449 CASE(dbcsr_conjugate_transpose) 00450 CALL dbcsr_init(matrix_left) 00451 CALL dbcsr_new_transposed (matrix_left, matrix_a,& 00452 shallow_data_copy=.FALSE., redistribute=.FALSE.,& 00453 transpose_distribution=.FALSE., error=error) 00454 CALL dbcsr_conjg(matrix_left, error=error) 00455 new_left = .TRUE. 00456 CASE DEFAULT 00457 CALL dbcsr_assert(.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, & 00458 routineN, "wrong transa_l = "//transa_l, __LINE__, error) 00459 END SELECT 00460 00461 SELECT CASE(transb_l) 00462 CASE(dbcsr_no_transpose) 00463 matrix_right = matrix_b 00464 new_right = .FALSE. 00465 CASE(dbcsr_transpose) 00466 CALL dbcsr_init(matrix_right) 00467 IF(dbcsr_get_matrix_type(matrix_b).EQ.dbcsr_type_antisymmetric) THEN 00468 ! 00469 ! For antisymmetric matrix, we need to do a hard copy 00470 ! shallow_data_copy=.TRUE. doesnt handle properly antisymm matrices 00471 CALL dbcsr_new_transposed (matrix_right, matrix_b,& 00472 shallow_data_copy=.FALSE., redistribute=.FALSE.,& 00473 transpose_distribution=.FALSE., error=error) 00474 ELSE 00475 CALL dbcsr_new_transposed (matrix_right, matrix_b,& 00476 shallow_data_copy=.TRUE., redistribute=.FALSE.,& 00477 transpose_distribution=.FALSE., error=error) 00478 ENDIF 00479 new_right = .TRUE. 00480 CASE(dbcsr_conjugate_transpose) 00481 CALL dbcsr_init(matrix_right) 00482 CALL dbcsr_new_transposed (matrix_right, matrix_b,& 00483 shallow_data_copy=.FALSE., redistribute=.FALSE.,& 00484 transpose_distribution=.FALSE., error=error) 00485 CALL dbcsr_conjg(matrix_right, error=error) 00486 new_right = .TRUE. 00487 CASE DEFAULT 00488 CALL dbcsr_assert(.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, & 00489 routineN, "wrong transb_l = "//transb_l, __LINE__, error) 00490 END SELECT 00491 00492 ! 00493 ! Ensure matrix compatibility. 00494 CALL dbcsr_assert (array_equality (dbcsr_row_block_offsets (matrix_c),& 00495 dbcsr_row_block_offsets (matrix_left)),& 00496 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00497 "C/A rows not equal", __LINE__, error=error) 00498 CALL dbcsr_assert (array_equality (dbcsr_col_block_offsets (matrix_c),& 00499 dbcsr_col_block_offsets (matrix_right)),& 00500 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00501 "C/B columns not equal", __LINE__, error=error) 00502 CALL dbcsr_assert (array_equality (dbcsr_col_block_offsets (matrix_left),& 00503 dbcsr_row_block_offsets (matrix_right)),& 00504 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00505 "A cols/B rows not equal", __LINE__, error=error) 00506 00507 ! 00508 ! No dense multiplication when filtering is used. 00509 use_dense_mult = .NOT. PRESENT (filter_eps) 00510 IF (mm_async) use_dense_mult = .FALSE. 00511 ! we skip dense multiply for (anti)symmetric matrices (slowdown for S/H * C) 00512 IF (use_dense_mult) THEN 00513 IF(dbcsr_has_symmetry (matrix_left) .OR. & 00514 dbcsr_has_symmetry(matrix_right)) THEN 00515 use_dense_mult = .FALSE. 00516 ELSE 00517 use_dense_mult = dbcsr_may_be_dense (matrix_left, make_dense_occ_thresh)& 00518 .AND. dbcsr_may_be_dense (matrix_right, make_dense_occ_thresh) 00519 ENDIF 00520 ENDIF 00521 ab_dense = use_dense_mult 00522 ! 00523 ! Submatrix selection 00524 f_row = 1 00525 l_row = dbcsr_nfullrows_total(matrix_c) 00526 f_col = 1 00527 l_col = dbcsr_nfullcols_total(matrix_c) 00528 f_k = 0 00529 l_k = 0 00530 IF (PRESENT (first_row)) THEN 00531 CALL dbcsr_assert(first_row .GE. 1& 00532 .AND. first_row .LE. dbcsr_nfullrows_total(matrix_c),& 00533 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00534 "Invalid first row specified", __LINE__, error) 00535 f_row = first_row 00536 ENDIF 00537 IF (PRESENT (last_row)) THEN 00538 CALL dbcsr_assert(last_row .GE. 1& 00539 .AND. last_row .LE. dbcsr_nfullrows_total(matrix_c)& 00540 .OR. last_row .LT. 1 ,& 00541 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00542 "Invalid last row specified", __LINE__, error) 00543 l_row = last_row 00544 ENDIF 00545 IF (PRESENT (first_column)) THEN 00546 CALL dbcsr_assert(first_column .GE. 1& 00547 .AND. first_column .LE. dbcsr_nfullcols_total(matrix_c),& 00548 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00549 "Invalid first col specified", __LINE__, error) 00550 f_col = first_column 00551 ENDIF 00552 IF (PRESENT (last_column)) THEN 00553 CALL dbcsr_assert(last_column .GE. 1& 00554 .AND. last_column .LE. dbcsr_nfullcols_total(matrix_c)& 00555 .OR. last_column .LT. 1,& 00556 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00557 "Invalid last column specified (C)", __LINE__, error) 00558 CALL dbcsr_assert(last_column .GE. 1& 00559 .AND. last_column .LE. dbcsr_nfullcols_total(matrix_right)& 00560 .OR. last_column .LT. 1,& 00561 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00562 "Invalid last column specified (B)", __LINE__, error) 00563 l_col = last_column 00564 ENDIF 00565 IF (PRESENT (first_k)) THEN 00566 CALL dbcsr_assert(first_k .GE. 1& 00567 .AND. first_k .LE. dbcsr_nfullcols_total(matrix_left),& 00568 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00569 "Invalid first k specified (A)", __LINE__, error) 00570 f_k = first_k 00571 ENDIF 00572 IF (PRESENT (last_k)) THEN 00573 CALL dbcsr_assert(last_k.GE. 1& 00574 .AND. last_k .LE. dbcsr_nfullcols_total(matrix_left)& 00575 .OR. last_k .LT. 1,& 00576 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 00577 "Invalid last k specified (A)", __LINE__, error) 00578 l_k = last_k 00579 ENDIF 00580 ! 00581 ! Now optimize the default submatrix selection values away 00582 IF (f_row .EQ. 1) f_row = 0 00583 IF (l_row .EQ. dbcsr_nfullrows_total(matrix_left)) l_row = 0 00584 IF (f_col .EQ. 1) f_col = 0 00585 ! The last column must be set if the right and product matrices 00586 ! differ. 00587 l_col = MIN (l_col, dbcsr_nfullcols_total(matrix_right)) 00588 l_col = MIN (l_col, dbcsr_nfullcols_total(matrix_c)) 00589 IF (f_col.LE.1.AND.& 00590 l_col .EQ. dbcsr_nfullcols_total(matrix_right) .AND. & 00591 dbcsr_nfullcols_total(matrix_right) .EQ.& 00592 dbcsr_nfullcols_total(matrix_c)) l_col = 0 00593 IF (f_k .EQ. 1) f_k = 0 00594 IF (l_k .EQ. dbcsr_nfullcols_total(matrix_left)) l_k = 0 00595 IF (.NOT. PRESENT(last_column) .AND.& 00596 .NOT. array_equality (dbcsr_col_block_sizes (matrix_right),& 00597 dbcsr_col_block_sizes (matrix_c))) THEN 00598 l_col = MIN (dbcsr_nfullcols_total(matrix_right),& 00599 dbcsr_nfullcols_total(matrix_c)) 00600 ENDIF 00601 CALL dbcsr_assert (f_row .LE. l_row, dbcsr_fatal_level,& 00602 dbcsr_wrong_args_error, routineN, "Last row smaller than first row", & 00603 __LINE__, error) 00604 CALL dbcsr_assert (f_col .LE. l_col, dbcsr_fatal_level,& 00605 dbcsr_wrong_args_error, routineN, "Last col smaller than first col", & 00606 __LINE__, error) 00607 00608 ! 00609 ! if we have limits we need to turn of make dense for the moment... 00610 !IF(ANY((/ f_row, l_row, f_col, l_col, f_k, l_k /).NE.0)) use_dense_mult = .FALSE. 00611 00612 ! 00613 ! 00614 ! Product data needs to be retained when 00615 ! * beta != 0; or 00616 ! * there is column limiting (l_col > 0) and the limiting column 00617 ! is less than the number of full columns in theproduct matrix 00618 keep_sparsity = .FALSE. 00619 IF (PRESENT (retain_sparsity)) keep_sparsity=retain_sparsity 00620 ! 00621 keep_product_data = keep_sparsity& 00622 .OR. .NOT. dbcsr_scalar_are_equal (beta, dbcsr_scalar_zero(beta%data_type))& 00623 .OR. (l_col .GT. 0 .AND. l_col .LT. dbcsr_nfullcols_total(matrix_c)) & 00624 .OR. (l_row .GT. 0 .AND. l_row .LT. dbcsr_nfullrows_total(matrix_c)) 00625 ! 00626 IF (.NOT. dbcsr_scalar_are_equal (beta, dbcsr_scalar_one(beta%data_type)) .AND. keep_product_data) THEN 00627 CALL dbcsr_scale (matrix_c, alpha_scalar=beta, & 00628 limits=(/f_row,l_row,f_col,l_col/), error=error) 00629 ENDIF 00630 ! 00631 ! The index of the product matrix is twiddled into canonical form 00632 ! if it is (anti)symmetric (i.e., rows and columns are where the 00633 ! row/column distributions say they are). Doing this in advance 00634 ! makes the local multiply more efficient. 00635 IF (dbcsr_has_symmetry (matrix_c)) THEN 00636 product_reindex = .TRUE. 00637 ELSE 00638 product_reindex = .FALSE. 00639 ENDIF 00640 ! Product can not be made dense; however, A & B may still be made 00641 ! dense unless previously determined otherwise. 00642 IF (product_reindex.OR.keep_sparsity) THEN 00643 use_dense_mult = .FALSE. 00644 ENDIF 00645 ! 00646 ! The thread distribution must reflect the current (possibly 00647 ! dense) distribution 00648 !CALL dbcsr_assert (dbcsr_distribution_has_threads(product_matrix%m%dist),& 00649 ! dbcsr_fatal_level, dbcsr_internal_error, routineN,& 00650 ! "Thread distribution not defined.", __LINE__, error=error) 00651 IF (.NOT. dbcsr_distribution_has_threads(matrix_c%m%dist)) THEN 00652 release_tdist = .TRUE. 00653 CALL dbcsr_distribution_make_threads (matrix_c%m%dist) 00654 ELSE 00655 release_tdist = .FALSE. 00656 ENDIF 00657 ! 00658 ! Create imaged distributions for the multiply. 00659 CALL dbcsr_create_image_dist (rdist_right, matrix_right%m%dist,& 00660 match_row_nbins = dbcsr_mp_npcols (dbcsr_distribution_mp (matrix_left%m%dist)),& 00661 match_col_nbins = dbcsr_mp_npcols (dbcsr_distribution_mp (matrix_c%m%dist)),& 00662 match_col_pdist = array_data (dbcsr_distribution_col_dist (matrix_c%m%dist))) 00663 CALL dbcsr_create_image_dist (rdist_left, matrix_left%m%dist,& 00664 match_row_pdist = array_data (dbcsr_distribution_row_dist (matrix_c%m%dist)),& 00665 match_row_nbins = dbcsr_mp_nprows (dbcsr_distribution_mp (matrix_c%m%dist)),& 00666 match_col_pdist = array_data (dbcsr_distribution_row_dist (rdist_right%i%main)),& 00667 match_col_idist = array_data (rdist_right%i%row_image),& 00668 match_col_nbins = dbcsr_mp_nprows (dbcsr_distribution_mp(matrix_right%m%dist))) 00669 IF (ab_dense) THEN 00670 CALL dbcsr_make_dists_dense (dbcsr_distribution (matrix_c),& 00671 rdist_left, rdist_right, dense_product_distribution,& 00672 dense_rdist_left, dense_rdist_right, .not.use_dense_mult,& 00673 m_map, k_vmap, n_map, matrix_c%m%row_blk_size, error=error) 00674 CALL make_sizes_dense (matrix_c%m%row_blk_size, m_map,& 00675 dbcsr_distribution_nrows (dense_product_distribution),& 00676 dense_row_sizes,& 00677 error=error) 00678 CALL make_sizes_dense (matrix_c%m%col_blk_size, n_map, & 00679 dbcsr_distribution_ncols (dense_product_distribution),& 00680 dense_col_sizes,& 00681 error=error) 00682 CALL make_sizes_dense (matrix_right%m%row_blk_size, k_vmap,& 00683 dbcsr_distribution_nrows (dense_rdist_right%i%main),& 00684 dense_k_sizes,& 00685 error=error) 00686 CALL dbcsr_init (dense_template_left) 00687 CALL dbcsr_create (dense_template_left, template=matrix_left,& 00688 dist=dense_rdist_left%i%main,& 00689 row_blk_size=dense_row_sizes, col_blk_size=dense_k_sizes,& 00690 error=error) 00691 CALL dbcsr_init (dense_template_right) 00692 CALL dbcsr_create (dense_template_right, template=matrix_right,& 00693 dist=dense_rdist_right%i%main,& 00694 row_blk_size=dense_k_sizes, col_blk_size=dense_col_sizes,& 00695 error=error) 00696 ENDIF 00697 ! 00698 CALL dbcsr_assert (use_dense_mult, "IMP", ab_dense,& 00699 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 00700 "Wrong logic when making dense matrices.", __LINE__, error=error) 00701 IF (use_dense_mult) THEN 00702 old_product_row_blk_offsets = matrix_c%m%row_blk_offset 00703 old_product_col_blk_offsets = matrix_c%m%col_blk_offset 00704 old_product_row_blk_sizes = matrix_c%m%row_blk_size 00705 old_product_col_blk_sizes = matrix_c%m%col_blk_size 00706 CALL array_hold (old_product_row_blk_offsets) 00707 CALL array_hold (old_product_col_blk_offsets) 00708 CALL array_hold (old_product_row_blk_sizes) 00709 CALL array_hold (old_product_col_blk_sizes) 00710 old_product_distribution = dbcsr_distribution (matrix_c) 00711 CALL dbcsr_distribution_hold (old_product_distribution) 00712 CALL dbcsr_init (product_matrix) 00713 CALL dbcsr_make_dense (matrix_c, product_matrix,& 00714 dense_product_distribution,& 00715 dense_row_sizes, dense_col_sizes,& 00716 m_map, n_map,& 00717 error=error) 00718 ELSE 00719 CALL dbcsr_init (product_matrix) 00720 CALL dbcsr_copy(product_matrix, matrix_c, shallow_data=.TRUE., error=error) 00721 ENDIF 00722 IF (ab_dense) THEN 00723 CALL dbcsr_distribution_release (dense_product_distribution) 00724 ENDIF 00725 ! 00726 ! if we need a special core gemm, we set it here 00727 ! 00728 ! initialize plasma 00729 IF (mm_driver .EQ. mm_driver_plasma) THEN 00730 CALL dbcsr_plasma_init(plasma_is_set, error) 00731 use_plasma = plasma_is_set.AND.use_dense_mult 00732 ELSE 00733 use_plasma = .FALSE. 00734 ENDIF 00735 ! 00736 IF (.TRUE. .OR. use_local_indexing) THEN 00737 ! This is needed to build the hash tables because they are 00738 ! locally indexed. 00739 CALL dbcsr_reset_locals (product_matrix, error=error) 00740 ENDIF 00741 IF (debug_mod) THEN 00742 WRITE(*,*)routineN//" Matrices ", dbcsr_get_matrix_type(matrix_a),& 00743 dbcsr_get_matrix_type(matrix_b), dbcsr_get_matrix_type(matrix_c) 00744 WRITE(*,*)routineN//" Matrices ", transa_l, transb_l, "keep", keep_product_data 00745 ENDIF 00746 IF (keep_product_data) THEN 00747 IF (product_reindex) THEN 00748 IF (debug_mod) WRITE(*,*)routineN//" Making canonical index" 00749 CALL dbcsr_make_index_canonical (product_matrix) 00750 ENDIF 00751 CALL dbcsr_assert (.NOT. ASSOCIATED (product_matrix%m%wms),& 00752 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 00753 "Product matrix should be finalized!", __LINE__, error=error) 00754 CALL dbcsr_make_untransposed_blocks (product_matrix, error=error) 00755 !$omp parallel if( .NOT. use_plasma ) & 00756 !$omp default (none) shared (product_matrix, error) 00757 ! For the multiply logic to work correctly, existing data must 00758 ! be added only after the index has been transformed into the 00759 ! canonical form. 00760 CALL dbcsr_add_wm_from_matrix(product_matrix, error=error) 00761 !$omp end parallel 00762 ELSE 00763 !$omp parallel if (.NOT. use_plasma) 00764 CALL dbcsr_work_create(product_matrix, work_mutable=.FALSE., error=error) 00765 !$omp end parallel 00766 ENDIF 00767 product_matrix%m%nze = 0 00768 product_matrix%m%row_p(:) = 0 00769 CALL dbcsr_data_set_size_referenced(product_matrix%m%data_area, 0) 00770 product_matrix%m%nblks = 0 00771 product_matrix%m%valid = .FALSE. 00772 ! 00773 NULLIFY (m2s_right) 00774 right_set_given = .FALSE. 00775 IF (PRESENT (right_set)) THEN 00776 IF (ASSOCIATED (right_set)) THEN 00777 right_set_given = .TRUE. 00778 m2s_right => right_set 00779 ENDIF 00780 ENDIF 00781 NULLIFY (m2s_left) 00782 left_set_given = .FALSE. 00783 IF (PRESENT (left_set)) THEN 00784 IF (ASSOCIATED (left_set)) THEN 00785 left_set_given = .TRUE. 00786 m2s_left => left_set 00787 ENDIF 00788 ENDIF 00789 ! 00790 ! Right images 00791 IF (.NOT. right_set_given) THEN 00792 ALLOCATE (m2s_right) 00793 IF (.NOT. dbcsr_scalar_are_equal (alpha, dbcsr_scalar_one(alpha%data_type))) THEN 00794 ! Copy and scale matrix B if alpha is not 1. 00795 CALL dbcsr_make_images (matrix_right, m2s_right, rdist_right,& 00796 predistribute="R", & 00797 data_memory_type = dbcsr_memory_default,& 00798 index_memory_type = dbcsr_memory_default,& 00799 no_copy_data=use_dense_mult, scale_value=alpha, error=error) 00800 ELSE 00801 CALL dbcsr_make_images (matrix_right, m2s_right, rdist_right,& 00802 predistribute="R", & 00803 data_memory_type = dbcsr_memory_default,& 00804 index_memory_type = dbcsr_memory_default,& 00805 no_copy_data=use_dense_mult, error=error) 00806 ENDIF 00807 ! Post-processing of images. 00808 DO i = 1, SIZE (m2s_right%mats,1) 00809 CALL dbcsr_reset_vlocals (m2s_right%mats(i,1), rdist_right, error=error) 00810 ! Crop if necessary 00811 IF (ANY ((/ f_k, l_k, f_col, l_col /) .NE. 0)) THEN 00812 CALL dbcsr_init (matrix_tmp) 00813 CALL dbcsr_crop_matrix (matrix_tmp, m2s_right%mats(i,1),& 00814 full_row_bounds=((/ f_k, l_k /)),& 00815 full_column_bounds=((/ f_col, l_col /)),& 00816 shallow_data = .FALSE., error=error) 00817 CALL dbcsr_release (m2s_right%mats(i,1)) 00818 CALL dbcsr_copy (m2s_right%mats(i,1), matrix_tmp, shallow_data=.TRUE.,& 00819 error=error) 00820 CALL dbcsr_release (matrix_tmp) 00821 CALL dbcsr_reset_vlocals (m2s_right%mats(i,1), rdist_right, error=error) 00822 ENDIF 00823 ENDDO 00824 IF (ab_dense) THEN 00825 CALL dbcsr_make_images_dense (m2s_right, dense_rdist_right, & 00826 row_map = k_vmap, col_map = n_map,& 00827 join_cols = use_dense_mult, join_rows=ab_dense, & 00828 new_template=dense_template_right, error=error) 00829 CALL dbcsr_image_dist_release (rdist_right, error=error) 00830 rdist_right = dense_rdist_right 00831 CALL dbcsr_image_dist_hold (rdist_right, error=error) 00832 DO i = 1, SIZE (m2s_right%mats,1) 00833 CALL dbcsr_reset_vlocals (m2s_right%mats(i,1), rdist_right, error=error) 00834 ENDDO 00835 ENDIF 00836 IF (use_local_indexing) THEN 00837 ! Convert to local-row index 00838 DO i = 1, SIZE (m2s_right%mats,1) 00839 CALL dbcsr_make_index_local_row(m2s_right%mats(i,1), error=error) 00840 ENDDO 00841 ENDIF 00842 IF (use_list_indexing) THEN 00843 ! Convert to list index 00844 DO i = 1, SIZE (m2s_right%mats,1) 00845 CALL dbcsr_make_index_list(m2s_right%mats(i,1), thread_redist=.FALSE.,& 00846 error=error) 00847 ENDDO 00848 ENDIF 00849 IF (use_local_indexing .AND. .NOT. use_list_indexing) THEN 00850 DO i = 1, SIZE (m2s_right%mats,1) 00851 CALL dbcsr_index_compact(m2s_right%mats(i,1), error=error) 00852 ENDDO 00853 ENDIF 00854 IF (PRESENT (right_set)) THEN 00855 right_set => m2s_right 00856 ENDIF 00857 ELSE 00858 m2s_right => right_set 00859 ENDIF 00860 IF (ab_dense) THEN 00861 CALL dbcsr_image_dist_release (dense_rdist_right, error=error) 00862 ENDIF 00863 ! 00864 ! Left images 00865 IF (.NOT. left_set_given) THEN 00866 ALLOCATE (m2s_left) 00867 CALL dbcsr_make_images (matrix_left, m2s_left, rdist_left,& 00868 predistribute="L", & 00869 data_memory_type = dbcsr_memory_default,& 00870 index_memory_type = dbcsr_memory_default,& 00871 no_copy_data=use_dense_mult, error=error) 00872 ! Post-processing of images. 00873 DO i = 1, SIZE (m2s_left%mats,2) 00874 CALL dbcsr_reset_vlocals (m2s_left%mats(1,i), rdist_left, error=error) 00875 ! Crop if necessary 00876 IF (ANY ((/ f_row, l_row, f_k, l_k /) .NE. 0)) THEN 00877 CALL dbcsr_init (matrix_tmp) 00878 CALL dbcsr_crop_matrix (matrix_tmp, m2s_left%mats(1,i),& 00879 full_row_bounds=((/ f_row, l_row /)),& 00880 full_column_bounds=((/ f_k, l_k /)),& 00881 shallow_data = .FALSE., error=error) 00882 CALL dbcsr_release (m2s_left%mats(1,i)) 00883 CALL dbcsr_copy (m2s_left%mats(1,i), matrix_tmp, shallow_data=.TRUE.,& 00884 error=error) 00885 CALL dbcsr_release (matrix_tmp) 00886 CALL dbcsr_reset_vlocals (m2s_left%mats(1,i), rdist_left, error=error) 00887 ENDIF 00888 ENDDO 00889 IF (ab_dense) THEN 00890 CALL dbcsr_make_images_dense (m2s_left, dense_rdist_left,& 00891 row_map = m_map, col_map = k_vmap,& 00892 join_rows = use_dense_mult, join_cols=ab_dense,& 00893 new_template=dense_template_left, error=error) 00894 CALL dbcsr_image_dist_release (rdist_left, error=error) 00895 rdist_left = dense_rdist_left 00896 CALL dbcsr_image_dist_hold (rdist_left, error=error) 00897 DO i = 1, SIZE (m2s_left%mats,2) 00898 CALL dbcsr_reset_vlocals (m2s_left%mats(1,i), rdist_left, error=error) 00899 ENDDO 00900 ENDIF 00901 IF (use_local_indexing) THEN 00902 ! Convert to local-row index 00903 DO i = 1, SIZE (m2s_left%mats,2) 00904 CALL dbcsr_make_index_local_row (m2s_left%mats(1,i), error=error) 00905 ENDDO 00906 END IF 00907 IF (use_list_indexing) THEN 00908 ! Convert to list index 00909 DO i = 1, SIZE (m2s_left%mats,2) 00910 CALL dbcsr_make_index_list (m2s_left%mats(1,i), thread_redist=.TRUE.,& 00911 error=error) 00912 ENDDO 00913 END IF 00914 IF (use_local_indexing .AND. .NOT. use_list_indexing) THEN 00915 DO i = 1, SIZE (m2s_left%mats,2) 00916 CALL dbcsr_index_compact (m2s_left%mats(1,i), error=error) 00917 ENDDO 00918 ENDIF 00919 IF (PRESENT (left_set)) THEN 00920 left_set => m2s_left 00921 ENDIF 00922 ELSE 00923 m2s_left => left_set 00924 ENDIF 00925 IF (ab_dense) THEN 00926 CALL dbcsr_image_dist_release (dense_rdist_left, error=error) 00927 ENDIF 00928 ! 00929 IF (ab_dense) THEN 00930 CALL array_release (k_vmap) 00931 CALL dbcsr_release (dense_template_left) 00932 CALL dbcsr_release (dense_template_right) 00933 CALL array_release (dense_row_sizes) 00934 CALL array_release (dense_col_sizes) 00935 CALL array_release (dense_k_sizes) 00936 ENDIF 00937 ! 00938 ! The limits were already used. Reset them. 00939 f_row = 0 ; l_row = 0 00940 f_col = 0 ; l_col = 0 00941 f_k = 0 ; l_k = 0 00942 ! 00943 ! Flush 00944 DO i = 1, SIZE (m2s_right%mats,1) 00945 CALL dbcsr_access_flush (m2s_right%mats(i,1), error=error) 00946 ENDDO 00947 DO i = 1, SIZE (m2s_left%mats,2) 00948 CALL dbcsr_access_flush (m2s_left%mats(1,i), error=error) 00949 ENDDO 00950 CALL dbcsr_access_flush (product_matrix, error=error) 00951 ! 00952 my_flop = 0 00953 CALL dbcsr_mm_cannon_multiply(m2s_left, m2s_right, product_matrix,& 00954 retain_sparsity=retain_sparsity,& 00955 filter_eps=filter_eps, error=error,& 00956 flop=my_flop) 00957 CALL dbcsr_finalize(product_matrix, error=error) 00958 IF (PRESENT (flop)) THEN 00959 ! return the average number of flops per MPI rank. Variance (which is fairly large) could be computed as well. 00960 comm = dbcsr_mp_group (dbcsr_distribution_mp (dbcsr_distribution (product_matrix))) 00961 numnodes = dbcsr_mp_numnodes (dbcsr_distribution_mp (dbcsr_distribution (product_matrix))) 00962 CALL mp_sum(my_flop,comm) 00963 flop = (my_flop + numnodes - 1) / numnodes 00964 ENDIF 00965 ! 00966 IF (new_left) CALL dbcsr_release (matrix_left) 00967 IF (new_right) CALL dbcsr_release (matrix_right) 00968 IF (release_tdist) THEN 00969 CALL dbcsr_distribution_no_threads (product_matrix%m%dist) 00970 ENDIF 00971 ! 00972 IF (.TRUE. .OR. use_local_indexing) & 00973 CALL dbcsr_release_locals (product_matrix, error=error) 00974 ! The index of the product matrix is reset to the CP2K form if it 00975 ! was previously set to the canonical form. 00976 IF (product_reindex) THEN 00977 IF (debug_mod) WRITE(*,*)routineN//" Making CP2K index" 00978 CALL dbcsr_make_index_canonical(product_matrix, cp2k=.TRUE.) 00979 ENDIF 00980 IF (use_dense_mult) THEN 00981 CALL dbcsr_release (matrix_c) 00982 CALL dbcsr_init (matrix_c) 00983 CALL dbcsr_make_undense(product_matrix, matrix_c,& 00984 old_product_distribution,& 00985 old_product_row_blk_offsets, old_product_col_blk_offsets,& 00986 old_product_row_blk_sizes, old_product_col_blk_sizes,& 00987 m_map, n_map, error=error) 00988 CALL dbcsr_release (product_matrix) 00989 CALL array_release (old_product_row_blk_offsets) 00990 CALL array_release (old_product_col_blk_offsets) 00991 CALL array_release (old_product_row_blk_sizes) 00992 CALL array_release (old_product_col_blk_sizes) 00993 CALL dbcsr_distribution_release (old_product_distribution) 00994 ELSE 00995 CALL dbcsr_release (matrix_c) 00996 CALL dbcsr_init (matrix_c) 00997 CALL dbcsr_copy (matrix_c, product_matrix, shallow_data=.TRUE., error=error) 00998 CALL dbcsr_release (product_matrix) 00999 ENDIF 01000 ! 01001 IF (.NOT. PRESENT (left_set)) THEN 01002 CALL dbcsr_destroy_array (m2s_left, error=error) 01003 DEALLOCATE (m2s_left) 01004 ENDIF 01005 CALL dbcsr_image_dist_release (rdist_left, error=error) 01006 IF (.NOT. PRESENT (right_set)) THEN 01007 CALL dbcsr_destroy_array (m2s_right, error=error) 01008 DEALLOCATE (m2s_right) 01009 ENDIF 01010 CALL dbcsr_image_dist_release (rdist_right, error=error) 01011 IF (ab_dense) THEN 01012 CALL array_release (m_map) 01013 CALL array_release (n_map) 01014 ENDIF 01015 ! 01016 ! if filtering is requested remove small blocks, unless the sparsity needs to be kept. 01017 ! 01018 IF (PRESENT (filter_eps) .AND. .NOT. keep_sparsity) THEN 01019 eps_any = dbcsr_scalar(filter_eps) 01020 CALL dbcsr_scalar_fill_all(eps_any) 01021 CALL dbcsr_scalar_set_type(eps_any, dbcsr_get_data_type(matrix_c)) 01022 CALL dbcsr_filter (matrix_c, eps_any, quick=.FALSE., error=error) 01023 ENDIF 01024 ! 01025 ! To support the canonical multiply (all non-transposed blocks), 01026 ! blocks may have to be transposed according to the CP2K 01027 ! triangular index. 01028 CALL dbcsr_make_untransposed_blocks (matrix_c, error=error) 01029 ! 01030 IF (debug_mod .OR. careful_mod) THEN 01031 IF (debug_mod) & 01032 WRITE(*,*)routineN//" Use dense mult, symm",& 01033 use_dense_mult, ab_dense, dbcsr_has_symmetry (matrix_c) 01034 CALL dbcsr_verify_matrix (matrix_c, error=error) 01035 IF (debug_mod) THEN 01036 cs = dbcsr_checksum (matrix_c, error=error) 01037 WRITE(*,*)routineN//" Product checksum", cs 01038 ENDIF 01039 ENDIF 01040 ! 01041 IF (mm_driver .EQ. mm_driver_plasma) THEN 01042 ! finalize special core gemm 01043 CALL dbcsr_plasma_finalize(error) 01044 ENDIF 01045 IF (.FALSE.) WRITE(*,*)"Finished with one multiplication." 01046 output_unit = default_output_unit 01047 CALL dbcsr_error_stop(error_handler, error) 01048 END SUBROUTINE dbcsr_multiply_anytype 01049 01050 01051 01052 SUBROUTINE dbcsr_multiply_s(transa, transb,& 01053 alpha, matrix_a, matrix_b, beta, matrix_c,& 01054 first_row, last_row, first_column, last_column, first_k, last_k,& 01055 retain_sparsity, left_set, right_set, filter_eps,& 01056 error, flop) 01057 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 01058 REAL(KIND=real_4), INTENT(IN) :: alpha 01059 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 01060 REAL(KIND=real_4), INTENT(IN) :: beta 01061 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 01062 INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, 01063 first_column, last_column, 01064 first_k, last_k 01065 LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity 01066 TYPE(dbcsr_2d_array_type), OPTIONAL, 01067 POINTER :: left_set, right_set 01068 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: filter_eps 01069 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01070 INTEGER(KIND=int_8), INTENT(OUT), 01071 OPTIONAL :: flop 01072 01073 CALL dbcsr_multiply_anytype(transa, transb,& 01074 dbcsr_scalar(alpha), matrix_a, matrix_b, dbcsr_scalar(beta), matrix_c,& 01075 first_row, last_row, first_column, last_column, first_k, last_k,& 01076 retain_sparsity, left_set=left_set, right_set=right_set,& 01077 filter_eps=filter_eps,& 01078 error=error, flop=flop) 01079 END SUBROUTINE dbcsr_multiply_s 01080 SUBROUTINE dbcsr_multiply_d(transa, transb,& 01081 alpha, matrix_a, matrix_b, beta, matrix_c,& 01082 first_row, last_row, first_column, last_column, first_k, last_k,& 01083 retain_sparsity, left_set, right_set, filter_eps,& 01084 error, flop) 01085 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 01086 REAL(KIND=real_8), INTENT(IN) :: alpha 01087 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 01088 REAL(KIND=real_8), INTENT(IN) :: beta 01089 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 01090 INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, 01091 first_column, last_column, 01092 first_k, last_k 01093 LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity 01094 TYPE(dbcsr_2d_array_type), OPTIONAL, 01095 POINTER :: left_set, right_set 01096 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: filter_eps 01097 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01098 INTEGER(KIND=int_8), INTENT(OUT), 01099 OPTIONAL :: flop 01100 01101 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_multiply_d', 01102 routineP = moduleN//':'//routineN 01103 01104 IF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND.& 01105 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4 .AND.& 01106 dbcsr_get_data_type(matrix_c) .EQ. dbcsr_type_real_4) THEN 01107 CALL dbcsr_multiply_anytype(transa, transb,& 01108 dbcsr_scalar(REAL(alpha,real_4)), matrix_a, matrix_b, 01109 dbcsr_scalar(REAL(beta,real_4)), matrix_c, 01110 first_row, last_row, first_column, last_column, first_k, last_k, 01111 retain_sparsity, left_set=left_set, right_set=right_set, 01112 filter_eps=filter_eps, 01113 error=error, flop=flop) 01114 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND.& 01115 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8 .AND.& 01116 dbcsr_get_data_type(matrix_c) .EQ. dbcsr_type_real_8) THEN 01117 CALL dbcsr_multiply_anytype(transa, transb,& 01118 dbcsr_scalar(alpha), matrix_a, matrix_b, dbcsr_scalar(beta), matrix_c,& 01119 first_row, last_row, first_column, last_column, first_k, last_k,& 01120 retain_sparsity, left_set=left_set, right_set=right_set,& 01121 filter_eps=filter_eps,& 01122 error=error, flop=flop) 01123 ELSE 01124 CALL dbcsr_assert (.FALSE., dbcsr_failure_level, dbcsr_internal_error,& 01125 routineP, "This combination of data types NYI",__LINE__, error) 01126 ENDIF 01127 END SUBROUTINE dbcsr_multiply_d 01128 SUBROUTINE dbcsr_multiply_c(transa, transb,& 01129 alpha, matrix_a, matrix_b, beta, matrix_c,& 01130 first_row, last_row, first_column, last_column, first_k, last_k,& 01131 retain_sparsity, left_set, right_set, filter_eps,& 01132 error, flop) 01133 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 01134 COMPLEX(KIND=real_4), INTENT(IN) :: alpha 01135 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 01136 COMPLEX(KIND=real_4), INTENT(IN) :: beta 01137 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 01138 INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, 01139 first_column, last_column, 01140 first_k, last_k 01141 LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity 01142 TYPE(dbcsr_2d_array_type), OPTIONAL, 01143 POINTER :: left_set, right_set 01144 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: filter_eps 01145 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01146 INTEGER(KIND=int_8), INTENT(OUT), 01147 OPTIONAL :: flop 01148 01149 CALL dbcsr_multiply_anytype(transa, transb,& 01150 dbcsr_scalar(alpha), matrix_a, matrix_b, dbcsr_scalar(beta), matrix_c,& 01151 first_row, last_row, first_column, last_column, first_k, last_k,& 01152 retain_sparsity, left_set=left_set, right_set=right_set,& 01153 filter_eps=filter_eps,& 01154 error=error, flop=flop) 01155 END SUBROUTINE dbcsr_multiply_c 01156 SUBROUTINE dbcsr_multiply_z(transa, transb,& 01157 alpha, matrix_a, matrix_b, beta, matrix_c,& 01158 first_row, last_row, first_column, last_column, first_k, last_k,& 01159 retain_sparsity, left_set, right_set, filter_eps,& 01160 error, flop) 01161 CHARACTER(LEN=1), INTENT(IN) :: transa, transb 01162 COMPLEX(KIND=real_8), INTENT(IN) :: alpha 01163 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 01164 COMPLEX(KIND=real_8), INTENT(IN) :: beta 01165 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 01166 INTEGER, INTENT(IN), OPTIONAL :: first_row, last_row, 01167 first_column, last_column, 01168 first_k, last_k 01169 LOGICAL, INTENT(IN), OPTIONAL :: retain_sparsity 01170 TYPE(dbcsr_2d_array_type), OPTIONAL, 01171 POINTER :: left_set, right_set 01172 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: filter_eps 01173 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01174 INTEGER(KIND=int_8), INTENT(OUT), 01175 OPTIONAL :: flop 01176 01177 CALL dbcsr_multiply_anytype(transa, transb,& 01178 dbcsr_scalar(alpha), matrix_a, matrix_b, dbcsr_scalar(beta), matrix_c,& 01179 first_row, last_row, first_column, last_column, first_k, last_k,& 01180 retain_sparsity, left_set=left_set, right_set=right_set,& 01181 filter_eps=filter_eps,& 01182 error=error, flop=flop) 01183 END SUBROUTINE dbcsr_multiply_z 01184 01185 ! ***************************************************************************** 01195 SUBROUTINE dbcsr_scale_anytype(matrix_a, alpha_scalar, limits, error) 01196 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01197 TYPE(dbcsr_scalar_type), INTENT(IN) :: alpha_scalar 01198 INTEGER, DIMENSION(4), INTENT(IN), 01199 OPTIONAL :: limits 01200 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01201 01202 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_anytype', 01203 routineP = moduleN//':'//routineN 01204 INTEGER, PARAMETER :: first_col_i = 3, 01205 first_row_i = 1, 01206 last_col_i = 4, last_row_i = 2 01207 01208 INTEGER :: a_col, a_col_size, a_row, a_row_size, col_offset, 01209 error_handler, row_offset, scale_col_offset, scale_col_size, 01210 scale_row_offset, scale_row_size 01211 INTEGER, DIMENSION(4) :: my_limits 01212 LOGICAL :: do_scale, has_limits, tr 01213 TYPE(dbcsr_data_obj) :: data_any 01214 TYPE(dbcsr_iterator) :: iter 01215 TYPE(dbcsr_scalar_type) :: one 01216 01217 ! --------------------------------------------------------------------------- 01218 01219 CALL dbcsr_error_set(routineN, error_handler, error) 01220 01221 ! Limits are only honored if the argument is present and any are 01222 ! non-zero. 01223 IF (PRESENT (limits)) THEN 01224 has_limits = ANY (limits(:) .NE. 0) 01225 ELSE 01226 has_limits = .FALSE. 01227 ENDIF 01228 my_limits(first_row_i) = 1 01229 my_limits(last_row_i) = dbcsr_nfullrows_total (matrix_a) 01230 my_limits(first_col_i) = 1 01231 my_limits(last_col_i) = dbcsr_nfullcols_total (matrix_a) 01232 IF (has_limits) THEN 01233 IF (limits(last_col_i) .NE. 0) THEN 01234 IF (debug_mod) THEN 01235 CALL dbcsr_assert (limits(last_col_i) .GE. 0, "AND",& 01236 limits(last_col_i) .LE. dbcsr_nfullcols_total (matrix_a),& 01237 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01238 "Specified last column is out of bounds.",__LINE__,error) 01239 ENDIF 01240 my_limits(last_col_i) = limits(last_col_i) 01241 ENDIF 01242 IF (limits(first_col_i) .NE. 0) THEN 01243 IF (debug_mod) THEN 01244 CALL dbcsr_assert (limits(first_col_i) .GE. 0, "AND",& 01245 limits(first_col_i) .LE. dbcsr_nfullcols_total (matrix_a),& 01246 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01247 "Specified first column is out of bounds.",__LINE__,error) 01248 ENDIF 01249 my_limits(first_col_i) = limits(first_col_i) 01250 ENDIF 01251 IF (limits(last_row_i) .NE. 0) THEN 01252 IF (debug_mod) THEN 01253 CALL dbcsr_assert (limits(last_row_i) .GE. 0, "AND",& 01254 limits(last_row_i) .LE. dbcsr_nfullrows_total (matrix_a),& 01255 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01256 "Specified last row is out of bounds.",__LINE__,error) 01257 ENDIF 01258 my_limits(last_row_i) = limits(last_row_i) 01259 ENDIF 01260 IF (limits(first_row_i) .NE. 0) THEN 01261 IF (debug_mod) THEN 01262 CALL dbcsr_assert (limits(first_row_i) .GE. 0, "AND",& 01263 limits(first_row_i) .LE. dbcsr_nfullrows_total (matrix_a),& 01264 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01265 "Specified first row is out of bounds.",__LINE__,error) 01266 ENDIF 01267 my_limits(first_row_i) = limits(first_row_i) 01268 ENDIF 01269 ENDIF 01270 ! 01271 ! quick return if possible 01272 one = dbcsr_scalar_one (dbcsr_scalar_get_type (alpha_scalar)) 01273 do_scale = .NOT. dbcsr_scalar_are_equal (alpha_scalar, one) 01274 ! 01275 ! let's go 01276 IF(do_scale) THEN 01277 !$OMP PARALLEL DEFAULT (none) & 01278 !$OMP PRIVATE (iter, data_any) & 01279 !$OMP PRIVATE (a_row, a_col, tr, a_row_size, a_col_size, & 01280 !$OMP row_offset, col_offset) & 01281 !$OMP PRIVATE (scale_row_size, scale_col_size,& 01282 !$OMP scale_row_offset, scale_col_offset) & 01283 !$OMP SHARED (matrix_a, my_limits, error, alpha_scalar) 01284 CALL dbcsr_data_init (data_any) 01285 CALL dbcsr_data_new (data_any, dbcsr_type_1d_to_2d (dbcsr_get_data_type (matrix_a))) 01286 CALL dbcsr_iterator_start(iter, matrix_a, read_only=.FALSE.,& 01287 contiguous_pointers = .FALSE., dynamic = .TRUE.,& 01288 dynamic_byrows = .TRUE., shared=.TRUE.) 01289 iterations: DO WHILE (dbcsr_iterator_blocks_left (iter)) 01290 CALL dbcsr_iterator_next_block(iter, a_row, a_col, data_any, tr,& 01291 row_size=a_row_size, col_size=a_col_size, & 01292 row_offset=row_offset, col_offset=col_offset) 01293 IF (a_row_size .GT. 0 .AND. a_col_size .GT. 0) THEN 01294 CALL frame_block_limit (a_row_size, row_offset,& 01295 my_limits(first_row_i), my_limits(last_row_i),& 01296 scale_row_size, scale_row_offset) 01297 CALL frame_block_limit (a_col_size, col_offset,& 01298 my_limits(first_col_i), my_limits(last_col_i),& 01299 scale_col_size, scale_col_offset) 01300 IF (tr) THEN 01301 CALL swap (scale_row_size, scale_col_size) 01302 CALL swap (scale_row_offset, scale_col_offset) 01303 ENDIF 01304 CALL dbcsr_block_scale (data_any, scale=alpha_scalar,& 01305 row_size=scale_row_size, col_size=scale_col_size,& 01306 lb=scale_row_offset, lb2=scale_col_offset, error=error) 01307 ENDIF 01308 ENDDO iterations 01309 CALL dbcsr_iterator_stop(iter) 01310 CALL dbcsr_data_clear_pointer (data_any) 01311 CALL dbcsr_data_release (data_any) 01312 !$OMP END PARALLEL 01313 ENDIF 01314 CALL dbcsr_error_stop(error_handler, error) 01315 END SUBROUTINE dbcsr_scale_anytype 01316 01317 01318 ! ***************************************************************************** 01328 ELEMENTAL SUBROUTINE frame_block_limit (block_size, block_offset,& 01329 first_limit, last_limit,& 01330 frame_size, frame_offset) 01331 INTEGER, INTENT(IN) :: block_size, block_offset, 01332 first_limit, last_limit 01333 INTEGER, INTENT(OUT) :: frame_size, frame_offset 01334 01335 INTEGER :: f, l 01336 01337 f = MAX (block_offset, first_limit) 01338 l = MIN (block_offset + block_size - 1, last_limit) 01339 frame_size = MAX(l - f + 1, 0) 01340 frame_offset = MIN(f - block_offset + 1, block_size) 01341 END SUBROUTINE frame_block_limit 01342 01343 01344 ! ***************************************************************************** 01352 SUBROUTINE dbcsr_scale_mat_any(matrix_a, alpha_matrix, side, error) 01353 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01354 TYPE(dbcsr_data_obj), INTENT(IN) :: alpha_matrix 01355 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: side 01356 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01357 01358 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_mat_any', 01359 routineP = moduleN//':'//routineN 01360 01361 INTEGER :: a_blk, a_col, a_col_size, 01362 a_nze, a_row, a_row_size, 01363 col_offset, error_handler, 01364 row_offset 01365 INTEGER, ALLOCATABLE, DIMENSION(:) :: m_offset 01366 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 01367 LOGICAL :: right, tr 01368 TYPE(dbcsr_data_obj) :: buffer, data_any, mat_tmp 01369 TYPE(dbcsr_iterator) :: iter 01370 01371 ! --------------------------------------------------------------------------- 01372 01373 CALL dbcsr_error_set(routineN, error_handler, error) 01374 01375 ! 01376 ! set vars 01377 right = .TRUE. 01378 IF(PRESENT(side)) THEN 01379 SELECT CASE(side) 01380 CASE('right');right = .TRUE. 01381 CASE('left' );right = .FALSE. 01382 CASE DEFAULT 01383 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level,& 01384 dbcsr_wrong_args_error, routineN,& 01385 "wrong side="//side,__LINE__,error) 01386 END SELECT 01387 ENDIF 01388 ! 01389 CALL dbcsr_data_init (buffer) 01390 ! Clean up these maxvalues 01391 CALL dbcsr_data_new (buffer, dbcsr_get_data_type(matrix_a),& 01392 data_size=MAXVAL(array_data (matrix_a%m%row_blk_size)) & 01393 *MAXVAL(array_data (matrix_a%m%col_blk_size))) 01394 CALL dbcsr_data_clear (buffer) ! Zero the buffer 01395 ! 01396 ! TODO: Clean up these offsets using the matrix-supplied ones 01397 row_blk_size => array_data (matrix_a%m%row_blk_size) 01398 col_blk_size => array_data (matrix_a%m%col_blk_size) 01399 ALLOCATE(m_offset(dbcsr_nblkrows_total(matrix_a))) 01400 m_offset(1) = 1 01401 IF(right) THEN 01402 DO a_col = 2,dbcsr_nblkcols_total(matrix_a) 01403 a_col_size = col_blk_size(a_col-1) 01404 m_offset(a_col) = m_offset(a_col-1) + a_col_size**2 01405 ENDDO 01406 ELSE 01407 DO a_row = 2,dbcsr_nblkrows_total(matrix_a) 01408 a_row_size = row_blk_size(a_row-1) 01409 m_offset(a_row) = m_offset(a_row-1) + a_row_size**2 01410 ENDDO 01411 ENDIF 01412 ! 01413 ! This data area holds the current column 01414 CALL dbcsr_data_init (mat_tmp) 01415 CALL dbcsr_data_new (mat_tmp, dbcsr_data_get_type (alpha_matrix)) 01416 ! 01417 ! let's go 01418 CALL dbcsr_data_init (data_any) 01419 CALL dbcsr_data_new (data_any, dbcsr_get_data_type (matrix_a)) 01420 CALL dbcsr_iterator_start(iter, matrix_a) 01421 iterations: DO WHILE (dbcsr_iterator_blocks_left (iter)) 01422 CALL dbcsr_iterator_next_block(iter, a_row, a_col, data_any, tr,& 01423 block_number=a_blk,& 01424 row_size=a_row_size, col_size=a_col_size, & 01425 row_offset=row_offset, col_offset=col_offset) 01426 a_nze = a_row_size * a_col_size 01427 IF (a_nze .EQ. 0) CYCLE ! Skip empty blocks 01428 ! 01429 ! let's scale 01430 IF(right) THEN 01431 !A = A * alpha 01432 mat_tmp = pointer_view (mat_tmp, alpha_matrix, m_offset(a_col)) 01433 CALL dbcsr_blas_gemm ('N','N',& 01434 a_row_size, a_col_size, a_col_size,& 01435 dbcsr_scalar_one (alpha_matrix%d%data_type),& 01436 data_any, a_row_size,& 01437 mat_tmp, a_col_size,& 01438 dbcsr_scalar_zero (alpha_matrix%d%data_type),& 01439 buffer, a_row_size) 01440 !CALL dgemm('N','N',a_row_size,a_col_size,a_col_size,& 01441 ! 1.0_dp,data_any%d%r_dp(1),a_row_size,& 01442 ! alpha_matrix%d%r_dp(m_offset(a_col)),a_col_size,& 01443 ! 0.0_dp,buffer%d%r_dp(1),a_row_size) 01444 CALL dbcsr_blas_copy (a_nze, buffer, 1, data_any, 1) 01445 ELSE 01446 !A = alpha * A 01447 mat_tmp = pointer_view (mat_tmp, alpha_matrix, m_offset(a_row)) 01448 CALL dbcsr_blas_gemm ('N','N',& 01449 a_row_size, a_col_size, a_row_size,& 01450 dbcsr_scalar_one (alpha_matrix%d%data_type),& 01451 mat_tmp, a_row_size,& 01452 data_any, a_row_size,& 01453 dbcsr_scalar_zero (alpha_matrix%d%data_type),& 01454 buffer, a_row_size) 01455 !CALL dgemm('N','N',a_row_size,a_col_size,a_row_size,& 01456 ! & 1.0_dp,alpha_matrix%d%r_dp(m_offset(a_col)),a_row_size,& 01457 ! & data_any%d%r_dp(1),a_row_size,& 01458 ! & 0.0_dp,buffer%d%r_dp(1),a_row_size) 01459 CALL dbcsr_blas_copy(a_nze, buffer, 1, data_any, 1) 01460 ENDIF 01461 ENDDO iterations 01462 CALL dbcsr_iterator_stop(iter) 01463 CALL dbcsr_data_clear_pointer (data_any) 01464 CALL dbcsr_data_release (data_any) 01465 DEALLOCATE(m_offset) 01466 !WRITE(*,*)"refs:", buffer%d%refcount, buffer%d%data_type 01467 CALL dbcsr_data_release (buffer) 01468 CALL dbcsr_data_clear_pointer (mat_tmp) 01469 CALL dbcsr_data_release (mat_tmp) 01470 CALL dbcsr_error_stop(error_handler, error) 01471 01472 END SUBROUTINE dbcsr_scale_mat_any 01473 01474 01475 ! ***************************************************************************** 01481 SUBROUTINE dbcsr_scale_by_vector_anytype(matrix_a, alpha, side, error) 01482 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01483 TYPE(dbcsr_data_obj), INTENT(IN), 01484 OPTIONAL :: alpha 01485 CHARACTER(LEN=*), INTENT(IN) :: side 01486 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01487 01488 CHARACTER(len=*), PARAMETER :: 01489 routineN = 'dbcsr_scale_by_vector_anytype', 01490 routineP = moduleN//':'//routineN 01491 01492 INTEGER :: a_blk, a_col, a_col_size, a_nze, a_row, a_row_size, 01493 col_offset, data_type, error_handler, i, row_offset 01494 LOGICAL :: right, tr 01495 TYPE(dbcsr_data_obj) :: data_any 01496 TYPE(dbcsr_iterator) :: iter 01497 01498 ! --------------------------------------------------------------------------- 01499 01500 CALL dbcsr_error_set(routineN, error_handler, error) 01501 01502 ! check that alpha and matrix have the same data type 01503 CALL dbcsr_assert (dbcsr_get_data_type (matrix_a).EQ.alpha%d%data_type, dbcsr_fatal_level,& 01504 dbcsr_wrong_args_error, routineN, "wrong data type matrix_a",__LINE__,error) 01505 ! 01506 ! set vars 01507 right = .TRUE. 01508 SELECT CASE(side) 01509 CASE('right');right = .TRUE. 01510 CASE('left' );right = .FALSE. 01511 CASE DEFAULT 01512 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, & 01513 routineN, "wrong side="//side,__LINE__,error) 01514 END SELECT 01515 ! 01516 ! let's go 01517 data_type = dbcsr_get_data_type (matrix_a) 01518 CALL dbcsr_data_init(data_any) 01519 CALL dbcsr_data_new(data_any, dbcsr_get_data_type (matrix_a)) 01520 CALL dbcsr_iterator_start(iter, matrix_a) 01521 DO WHILE (dbcsr_iterator_blocks_left (iter)) 01522 CALL dbcsr_iterator_next_block(iter, a_row, a_col, data_any, tr,& 01523 block_number=a_blk,& 01524 row_size=a_row_size, col_size=a_col_size, & 01525 row_offset=row_offset, col_offset=col_offset) 01526 a_nze = a_row_size * a_col_size 01527 IF (a_nze .EQ. 0) CYCLE ! Skip empty blocks 01528 ! 01529 ! let's scale 01530 IF(right) THEN 01531 DO i = 1,a_col_size 01532 SELECT CASE (data_type) 01533 CASE (dbcsr_type_real_4) 01534 CALL sscal(a_row_size,alpha%d%r_sp(col_offset+i-1),& 01535 data_any%d%r_sp((i-1)*a_row_size+1),1) 01536 CASE (dbcsr_type_real_8) 01537 CALL dscal(a_row_size,alpha%d%r_dp(col_offset+i-1),& 01538 data_any%d%r_dp((i-1)*a_row_size+1),1) 01539 CASE (dbcsr_type_complex_4) 01540 CALL cscal(a_row_size,alpha%d%c_sp(col_offset+i-1),& 01541 data_any%d%c_sp((i-1)*a_row_size+1),1) 01542 CASE (dbcsr_type_complex_8) 01543 CALL zscal(a_row_size,alpha%d%c_dp(col_offset+i-1),& 01544 data_any%d%c_dp((i-1)*a_row_size+1),1) 01545 END SELECT 01546 ENDDO 01547 ELSE 01548 DO i = 1,a_row_size 01549 SELECT CASE (data_type) 01550 CASE (dbcsr_type_real_4) 01551 CALL sscal(a_col_size,alpha%d%r_sp(row_offset+i-1),& 01552 data_any%d%r_sp(i),a_col_size) 01553 CASE (dbcsr_type_real_8) 01554 CALL dscal(a_col_size,alpha%d%r_dp(row_offset+i-1),& 01555 data_any%d%r_dp(i),a_col_size) 01556 CASE (dbcsr_type_complex_4) 01557 CALL cscal(a_col_size,alpha%d%c_sp(row_offset+i-1),& 01558 data_any%d%c_sp(i),a_col_size) 01559 CASE (dbcsr_type_complex_8) 01560 CALL zscal(a_col_size,alpha%d%c_dp(row_offset+i-1),& 01561 data_any%d%c_dp(i),a_col_size) 01562 END SELECT 01563 ENDDO 01564 ENDIF 01565 ENDDO 01566 CALL dbcsr_iterator_stop(iter) 01567 CALL dbcsr_data_clear_pointer (data_any) 01568 CALL dbcsr_data_release (data_any) 01569 CALL dbcsr_error_stop(error_handler, error) 01570 01571 END SUBROUTINE dbcsr_scale_by_vector_anytype 01572 01573 01574 ! ***************************************************************************** 01580 SUBROUTINE dbcsr_set_anytype(matrix, alpha, error) 01581 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 01582 TYPE(dbcsr_scalar_type), INTENT(IN) :: alpha 01583 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01584 01585 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_anytype', 01586 routineP = moduleN//':'//routineN 01587 01588 INTEGER :: blk, col, error_handler, row 01589 LOGICAL :: tr 01590 TYPE(dbcsr_data_obj) :: data_block 01591 TYPE(dbcsr_iterator) :: iter 01592 01593 ! --------------------------------------------------------------------------- 01594 ! 01595 01596 CALL dbcsr_error_set(routineN, error_handler, error) 01597 CALL dbcsr_data_init (data_block) 01598 CALL dbcsr_data_new (data_block, dbcsr_get_data_type (matrix)) 01599 CALL dbcsr_iterator_start(iter, matrix) 01600 DO WHILE (dbcsr_iterator_blocks_left(iter)) 01601 CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr, blk) 01602 CALL dbcsr_data_clear (data_block, value=alpha) 01603 ENDDO 01604 CALL dbcsr_iterator_stop(iter) 01605 CALL dbcsr_data_clear_pointer (data_block) 01606 CALL dbcsr_data_release (data_block) 01607 CALL dbcsr_error_stop(error_handler, error) 01608 ! 01609 END SUBROUTINE dbcsr_set_anytype 01610 01611 ! ***************************************************************************** 01616 SUBROUTINE dbcsr_conjg(matrix, error) 01617 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 01618 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01619 01620 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_conjg', 01621 routineP = moduleN//':'//routineN 01622 01623 INTEGER :: blk, col, data_type, 01624 error_handler, row 01625 LOGICAL :: tr 01626 TYPE(dbcsr_data_obj) :: data_any 01627 TYPE(dbcsr_iterator) :: iter 01628 01629 ! --------------------------------------------------------------------------- 01630 ! 01631 01632 CALL dbcsr_error_set(routineN, error_handler, error) 01633 data_type = dbcsr_get_data_type(matrix) 01634 CALL dbcsr_data_init (data_any) 01635 CALL dbcsr_data_new (data_any, data_type) 01636 CALL dbcsr_iterator_start(iter, matrix) 01637 DO WHILE (dbcsr_iterator_blocks_left(iter)) 01638 CALL dbcsr_iterator_next_block(iter, row, col, data_any, tr, blk) 01639 SELECT CASE (data_type) 01640 CASE (dbcsr_type_complex_4) 01641 data_any%d%c_sp = CONJG(data_any%d%c_sp) 01642 CASE (dbcsr_type_complex_8) 01643 data_any%d%c_dp = CONJG(data_any%d%c_dp) 01644 CASE DEFAULT 01645 ! needed for g95 01646 END SELECT 01647 ENDDO 01648 CALL dbcsr_iterator_stop(iter) 01649 CALL dbcsr_data_clear_pointer(data_any) 01650 CALL dbcsr_data_release(data_any) 01651 CALL dbcsr_error_stop(error_handler, error) 01652 END SUBROUTINE dbcsr_conjg 01653 01654 ! ***************************************************************************** 01663 SUBROUTINE dbcsr_add_anytype(matrix_a, matrix_b, alpha_scalar, beta_scalar, error) 01664 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01665 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01666 TYPE(dbcsr_scalar_type), INTENT(IN), 01667 OPTIONAL :: alpha_scalar, beta_scalar 01668 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01669 01670 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_anytype', 01671 routineP = moduleN//':'//routineN 01672 01673 INTEGER :: blk, col, data_type_a, 01674 data_type_b, error_handler, 01675 row, size_a, size_b 01676 LOGICAL :: do_scale, tr 01677 TYPE(dbcsr_data_obj) :: data_block 01678 TYPE(dbcsr_iterator) :: iter 01679 TYPE(dbcsr_scalar_type) :: my_alpha_scalar, 01680 my_beta_scalar 01681 01682 ! --------------------------------------------------------------------------- 01683 01684 CALL dbcsr_error_set(routineN, error_handler, error) 01685 CALL dbcsr_assert (dbcsr_valid_index (matrix_a), dbcsr_fatal_level,& 01686 dbcsr_caller_error, routineN, "Invalid matrix", __LINE__, error=error) 01687 01688 IF((dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_symmetric.OR.& 01689 dbcsr_get_matrix_type(matrix_a).EQ.dbcsr_type_antisymmetric).NEQV. & 01690 (dbcsr_get_matrix_type(matrix_b).EQ.dbcsr_type_symmetric.OR.& 01691 dbcsr_get_matrix_type(matrix_b).EQ.dbcsr_type_antisymmetric)) THEN 01692 CALL dbcsr_assert (.FALSE.,dbcsr_fatal_level,& 01693 dbcsr_unimplemented_error_nr, routineN, "Summing general with symmetric matrix NYI",& 01694 __LINE__,error) 01695 ENDIF 01696 01697 data_type_a = dbcsr_get_data_type(matrix_a) 01698 data_type_b = dbcsr_get_data_type(matrix_b) 01699 ! 01700 my_alpha_scalar = dbcsr_scalar_one (data_type_a) 01701 IF(PRESENT(alpha_scalar)) my_alpha_scalar = alpha_scalar 01702 my_beta_scalar = dbcsr_scalar_one (data_type_b) 01703 IF(PRESENT(beta_scalar)) my_beta_scalar = beta_scalar 01704 ! 01705 ! let's go 01706 CALL dbcsr_assert (dbcsr_nblkrows_total(matrix_a).EQ.dbcsr_nblkrows_total(matrix_b), & 01707 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, "matrices not consistent",__LINE__,error) 01708 01709 do_scale = dbcsr_scalar_are_equal (& 01710 my_beta_scalar, dbcsr_scalar_one (data_type_b)) 01711 01712 CALL dbcsr_scale(matrix_a, alpha_scalar=my_alpha_scalar, error=error) 01713 01714 ! Pre-size work arrays of matrix_a to avoid continuous reallocation. 01715 size_a = dbcsr_get_data_size_referenced (matrix_a) 01716 size_b = dbcsr_get_data_size_referenced (matrix_b) 01717 IF(.NOT.dbcsr_scalar_are_equal(my_beta_scalar,& 01718 dbcsr_scalar_zero(data_type_b))) THEN 01719 !$OMP PARALLEL DEFAULT (none) & 01720 !$OMP PRIVATE (iter, data_block) & 01721 !$OMP PRIVATE (row, col, tr, blk) & 01722 !$OMP SHARED (matrix_a, matrix_b, data_type_b, size_b, size_a) & 01723 !$OMP SHARED (do_scale, my_beta_scalar) & 01724 !$OMP SHARED (error) 01725 IF (size_b .GT. size_a .AND. matrix_b%m%nblks .GT. matrix_a%m%nblks) THEN 01726 CALL dbcsr_work_create (matrix_a,& 01727 nblks_guess = matrix_b%m%nblks - matrix_a%m%nblks,& 01728 sizedata_guess = size_b - size_a,& 01729 work_mutable = .FALSE., error=error) 01730 ELSE 01731 CALL dbcsr_work_create (matrix_a,& 01732 work_mutable = .FALSE., error=error) 01733 ENDIF 01734 !$OMP BARRIER 01735 CALL dbcsr_data_init (data_block) 01736 CALL dbcsr_data_new (data_block, data_type_b) 01737 CALL dbcsr_iterator_start(iter, matrix_b,& 01738 shared = .TRUE., read_only = .TRUE., contiguous_pointers=.FALSE.,& 01739 dynamic = .TRUE., dynamic_byrows = .TRUE.) 01740 01741 DO WHILE (dbcsr_iterator_blocks_left(iter)) 01742 01743 CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr, blk) 01744 01745 IF (do_scale) THEN 01746 CALL dbcsr_put_block(matrix_a, row, col, data_block, tr,& 01747 summation=.TRUE.) 01748 ELSE 01749 CALL dbcsr_put_block(matrix_a, row, col, data_block, tr,& 01750 summation=.TRUE., scale=my_beta_scalar) 01751 ENDIF 01752 01753 ENDDO 01754 01755 CALL dbcsr_iterator_stop(iter) 01756 CALL dbcsr_finalize (matrix_a, error=error) 01757 CALL dbcsr_data_clear_pointer (data_block) 01758 CALL dbcsr_data_release (data_block) 01759 !$OMP END PARALLEL 01760 ENDIF 01761 CALL dbcsr_error_stop(error_handler, error) 01762 END SUBROUTINE dbcsr_add_anytype 01763 01764 ! ***************************************************************************** 01777 SUBROUTINE dbcsr_add_reserved(matrix_a, matrix_b,& 01778 alpha_scalar, beta_scalar, error) 01779 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01780 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01781 TYPE(dbcsr_scalar_type), INTENT(IN), 01782 OPTIONAL :: alpha_scalar, beta_scalar 01783 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01784 01785 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_reserved', 01786 routineP = moduleN//':'//routineN 01787 01788 INTEGER :: blk, col, data_type_a, 01789 data_type_b, error_handler, 01790 nblkrows, nblks, row 01791 INTEGER, ALLOCATABLE, DIMENSION(:) :: b_row_i 01792 LOGICAL :: do_scale, tr 01793 TYPE(dbcsr_data_obj) :: data_block 01794 TYPE(dbcsr_iterator) :: iter 01795 TYPE(dbcsr_scalar_type) :: my_alpha_scalar, 01796 my_beta_scalar 01797 01798 ! --------------------------------------------------------------------------- 01799 01800 CALL dbcsr_error_set(routineN, error_handler, error) 01801 ! Checks for validity 01802 CALL dbcsr_assert (dbcsr_valid_index (matrix_a),& 01803 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01804 "Target matrix A must be valid.", __LINE__, error) 01805 CALL dbcsr_assert (dbcsr_valid_index (matrix_b),& 01806 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01807 "Source matrix B must be valid.", __LINE__, error) 01808 ! Shortcuts 01809 data_type_a = dbcsr_get_data_type(matrix_a) 01810 data_type_b = dbcsr_get_data_type(matrix_b) 01811 ! Process arguments 01812 my_alpha_scalar = dbcsr_scalar_one (data_type_a) 01813 IF(PRESENT(alpha_scalar)) my_alpha_scalar = alpha_scalar 01814 my_beta_scalar = dbcsr_scalar_one (data_type_b) 01815 IF(PRESENT(beta_scalar)) my_beta_scalar = beta_scalar 01816 ! Scale target matrix 01817 do_scale = dbcsr_scalar_are_equal (& 01818 my_beta_scalar, dbcsr_scalar_one (data_type_b)) 01819 CALL dbcsr_scale(matrix_a, alpha_scalar=my_alpha_scalar, error=error) 01820 ! 01821 IF(.NOT.dbcsr_scalar_are_equal(my_beta_scalar,& 01822 dbcsr_scalar_zero(data_type_b))) THEN 01823 ! Reserve the blocks to be added 01824 nblks = dbcsr_get_num_blocks (matrix_b) 01825 nblkrows = dbcsr_nblkrows_total (matrix_b) 01826 ALLOCATE (b_row_i(nblks)) 01827 CALL dbcsr_expand_row_index (matrix_b%m%row_p, b_row_i, nblkrows, nblks) 01828 CALL dbcsr_reserve_blocks (matrix_a, b_row_i, matrix_b%m%col_i,& 01829 error=error) 01830 DEALLOCATE (b_row_i) 01831 ! Now add the blocks 01832 CALL dbcsr_data_init (data_block) 01833 CALL dbcsr_data_new (data_block, data_type_b) 01834 CALL dbcsr_iterator_start(iter, matrix_b) 01835 DO WHILE (dbcsr_iterator_blocks_left(iter)) 01836 CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr, blk) 01837 IF (do_scale) THEN 01838 CALL dbcsr_put_block(matrix_a, row, col, data_block, tr,& 01839 summation=.TRUE.) 01840 ELSE 01841 CALL dbcsr_put_block(matrix_a, row, col, data_block, tr,& 01842 summation=.TRUE., scale=my_beta_scalar) 01843 ENDIF 01844 ENDDO 01845 CALL dbcsr_assert (dbcsr_valid_index (matrix_a),& 01846 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 01847 "Block reservations seem to be incompletely done.", __LINE__,& 01848 error=error) 01849 CALL dbcsr_iterator_stop(iter) 01850 CALL dbcsr_data_clear_pointer (data_block) 01851 CALL dbcsr_data_release (data_block) 01852 !CALL dbcsr_finalize (matrix_a, error=error) 01853 ENDIF 01854 CALL dbcsr_error_stop(error_handler, error) 01855 END SUBROUTINE dbcsr_add_reserved 01856 01857 01859 SUBROUTINE dbcsr_add_d(matrix_a, matrix_b, alpha_scalar, beta_scalar, error) 01860 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01861 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01862 REAL(real_8), INTENT(IN) :: alpha_scalar, beta_scalar 01863 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01864 01865 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_d', 01866 routineP = moduleN//':'//routineN 01867 01868 INTEGER :: error_handler 01869 01870 CALL dbcsr_error_set(routineN, error_handler, error) 01871 IF( dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & 01872 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8) THEN 01873 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01874 alpha_scalar=dbcsr_scalar(alpha_scalar),& 01875 beta_scalar=dbcsr_scalar(beta_scalar), error=error) 01876 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & 01877 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN 01878 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01879 alpha_scalar=dbcsr_scalar(REAL(alpha_scalar,real_4)), 01880 beta_scalar=dbcsr_scalar(REAL(beta_scalar,real_4)), error=error) 01881 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & 01882 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN 01883 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01884 alpha_scalar=dbcsr_scalar(alpha_scalar),& 01885 beta_scalar=dbcsr_scalar(REAL(beta_scalar,real_4)), error=error) 01886 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & 01887 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8) THEN 01888 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01889 alpha_scalar=dbcsr_scalar(REAL(alpha_scalar,real_4)), 01890 beta_scalar=dbcsr_scalar(beta_scalar), error=error) 01891 ELSE 01892 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01893 "Invalid combination of data type, NYI",__LINE__,error) 01894 ENDIF 01895 CALL dbcsr_error_stop(error_handler, error) 01896 END SUBROUTINE dbcsr_add_d 01897 01898 SUBROUTINE dbcsr_add_s(matrix_a, matrix_b, alpha_scalar, beta_scalar, error) 01899 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01900 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01901 REAL(real_4), INTENT(IN) :: alpha_scalar, beta_scalar 01902 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01903 01904 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_s', 01905 routineP = moduleN//':'//routineN 01906 01907 INTEGER :: error_handler 01908 01909 CALL dbcsr_error_set(routineN, error_handler, error) 01910 IF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & 01911 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN 01912 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01913 alpha_scalar=dbcsr_scalar(alpha_scalar),& 01914 beta_scalar=dbcsr_scalar(beta_scalar), error=error) 01915 ELSE 01916 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01917 "Invalid combination of data type, NYI",__LINE__,error) 01918 ENDIF 01919 CALL dbcsr_error_stop(error_handler, error) 01920 END SUBROUTINE dbcsr_add_s 01921 01922 SUBROUTINE dbcsr_add_z(matrix_a, matrix_b, alpha_scalar, beta_scalar, error) 01923 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01924 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01925 COMPLEX(real_8), INTENT(IN) :: alpha_scalar, beta_scalar 01926 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01927 01928 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_z', 01929 routineP = moduleN//':'//routineN 01930 01931 INTEGER :: error_handler 01932 01933 CALL dbcsr_error_set(routineN, error_handler, error) 01934 IF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_8 .AND. & 01935 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_8) THEN 01936 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01937 alpha_scalar=dbcsr_scalar(alpha_scalar),& 01938 beta_scalar=dbcsr_scalar(beta_scalar), error=error) 01939 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_4 .AND. & 01940 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_4) THEN 01941 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01942 alpha_scalar=dbcsr_scalar(CMPLX(alpha_scalar,KIND=real_4)),& 01943 beta_scalar=dbcsr_scalar(CMPLX(beta_scalar,KIND=real_4)), error=error) 01944 ELSE 01945 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01946 "Invalid combination of data type, NYI",__LINE__,error) 01947 ENDIF 01948 CALL dbcsr_error_stop(error_handler, error) 01949 END SUBROUTINE dbcsr_add_z 01950 01951 SUBROUTINE dbcsr_add_c(matrix_a, matrix_b, alpha_scalar, beta_scalar, error) 01952 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 01953 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 01954 COMPLEX(real_4), INTENT(IN) :: alpha_scalar, beta_scalar 01955 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 01956 01957 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_c', 01958 routineP = moduleN//':'//routineN 01959 01960 INTEGER :: error_handler 01961 01962 CALL dbcsr_error_set(routineN, error_handler, error) 01963 IF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_4 .AND. & 01964 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_4) THEN 01965 CALL dbcsr_add_anytype(matrix_a, matrix_b,& 01966 alpha_scalar=dbcsr_scalar(alpha_scalar),& 01967 beta_scalar=dbcsr_scalar(beta_scalar), error=error) 01968 ELSE 01969 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 01970 "Invalid combination of data type, NYI",__LINE__,error) 01971 ENDIF 01972 CALL dbcsr_error_stop(error_handler, error) 01973 END SUBROUTINE dbcsr_add_c 01974 01975 01976 ! ***************************************************************************** 01995 SUBROUTINE dbcsr_reduce (matrix, reduced,& 01996 reduction_target, reduce_rows, reduce_columns, & 01997 operation, error) 01998 TYPE(dbcsr_obj), INTENT(IN), TARGET :: matrix 01999 TYPE(dbcsr_obj), INTENT(INOUT), TARGET :: reduced 02000 INTEGER, INTENT(IN) :: reduction_target 02001 LOGICAL, INTENT(IN) :: reduce_rows, reduce_columns 02002 CHARACTER(LEN=*), INTENT(IN) :: operation 02003 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02004 02005 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reduce', 02006 routineP = moduleN//':'//routineN 02007 INTEGER, PARAMETER :: i_d = 2, i_i = 1 02008 REAL(kind=dp), PARAMETER :: comp_eps = 1.0E-06_dp 02009 02010 CHARACTER :: red_type 02011 INTEGER :: data_type, dst_p, error_handler, half_spread, most_spread, 02012 mp_group, my_norm_node, mynode, nsteps, numnodes, spread, src_p, step, 02013 tag 02014 INTEGER, DIMENSION(2) :: get_sizes, my_sizes 02015 LOGICAL :: i_recv, i_send 02016 REAL(kind=dp) :: cs1, cs2 02017 TYPE(dbcsr_data_obj) :: data_ptr 02018 TYPE(dbcsr_mp_obj) :: mp_obj 02019 TYPE(dbcsr_obj), POINTER :: mat_comm 02020 TYPE(dbcsr_obj), TARGET :: mat_recv 02021 02022 ! --------------------------------------------------------------------------- 02023 02024 CALL dbcsr_error_set (routineN, error_handler, error) 02025 CALL dbcsr_assert (dbcsr_valid_index (matrix%m), dbcsr_fatal_level,& 02026 dbcsr_caller_error, routineN, "Source matrix not initialized.",& 02027 __LINE__, error=error) 02028 CALL dbcsr_assert (operation, "EQ", "+", dbcsr_fatal_level,& 02029 dbcsr_unimplemented_error_nr, routineN,& 02030 "Reduction operation not yet implemented.", __LINE__, error=error) 02031 CALL dbcsr_access_flush (matrix, error=error) 02032 IF (reduce_rows .AND. reduce_columns) THEN 02033 red_type = dbcsr_repl_full 02034 ELSEIF (reduce_rows .AND. .NOT. reduce_columns) THEN 02035 red_type = dbcsr_repl_row 02036 ELSEIF (reduce_columns .AND. .NOT. reduce_rows) THEN 02037 red_type = dbcsr_repl_col 02038 ELSE 02039 red_type = dbcsr_repl_none 02040 CALL dbcsr_assert(reduce_rows, "OR", reduce_columns, & 02041 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02042 "At least one of row or column reduction must be specified",& 02043 __LINE__, error=error) 02044 ENDIF 02045 ! Remember stuff for easier handling. 02046 mp_obj = dbcsr_distribution_mp (dbcsr_distribution (matrix)) 02047 data_type = dbcsr_get_data_type (matrix) 02048 ! Set communication variables. 02049 SELECT CASE (red_type) 02050 CASE (dbcsr_repl_full) 02051 numnodes = dbcsr_mp_numnodes (mp_obj) 02052 mp_group = dbcsr_mp_group (mp_obj) 02053 mynode = dbcsr_mp_mynode (mp_obj) 02054 CASE (dbcsr_repl_row) 02055 numnodes = dbcsr_mp_npcols (mp_obj) 02056 CALL dbcsr_mp_grid_setup (mp_obj, force=.TRUE.) 02057 mp_group = dbcsr_mp_my_row_group (mp_obj) 02058 mynode = dbcsr_mp_mypcol (mp_obj) 02059 CASE (dbcsr_repl_col) 02060 numnodes = dbcsr_mp_nprows (mp_obj) 02061 CALL dbcsr_mp_grid_setup (mp_obj, force=.TRUE.) 02062 mp_group = dbcsr_mp_my_col_group (mp_obj) 02063 mynode = dbcsr_mp_myprow (mp_obj) 02064 CASE (dbcsr_repl_none) 02065 numnodes = 1 02066 mp_group = dbcsr_mp_group (mp_obj) 02067 mynode = dbcsr_mp_mynode (mp_obj) 02068 END SELECT 02069 ! Check that we actually have the row/column communicators when 02070 ! doing row or column-limited reduction. 02071 CALL dbcsr_assert (red_type.EQ.dbcsr_repl_row & 02072 .OR. red_type.EQ.dbcsr_repl_col,& 02073 "IMP", dbcsr_mp_has_subgroups (mp_obj), dbcsr_fatal_level,& 02074 dbcsr_unimplemented_error_nr, routineN,& 02075 "Only full reduction supported when subcommunicators are turned off.",& 02076 __LINE__, error=error) 02077 ! Check the reduction_target is sane 02078 CALL dbcsr_assert (reduction_target, "GE", 0, dbcsr_fatal_level,& 02079 dbcsr_wrong_args_error, routineN, "Can not have non-negative target.",& 02080 __LINE__, error=error) 02081 CALL dbcsr_assert (reduction_target, "LT", numnodes, dbcsr_fatal_level,& 02082 dbcsr_wrong_args_error, routineN, "Target must fit process grid.", & 02083 __LINE__, error=error) 02084 ! 02085 ! Send/receive buffer setup 02086 CALL dbcsr_init (mat_recv) 02087 CALL dbcsr_create (reduced, "Reduced "//matrix%m%name,& 02088 dbcsr_distribution (matrix), dbcsr_get_matrix_type (matrix),& 02089 dbcsr_row_block_sizes (matrix), dbcsr_col_block_sizes (matrix),& 02090 data_type=data_type, error=error) 02091 CALL dbcsr_finalize (reduced, error=error) 02092 CALL dbcsr_create (mat_recv, "Recv buffer for "//matrix%m%name,& 02093 dbcsr_distribution (matrix), dbcsr_get_matrix_type (matrix),& 02094 dbcsr_row_block_sizes (matrix), dbcsr_col_block_sizes (matrix),& 02095 data_type=data_type, data_memory_type=dbcsr_memory_MPI,& 02096 index_memory_type=dbcsr_memory_MPI, error=error) 02097 CALL dbcsr_data_init (data_ptr) 02098 CALL dbcsr_data_new (data_ptr, data_type) 02099 ! 02100 IF (debug_mod) THEN 02101 cs1 = dbcsr_checksum (matrix, error=error) 02102 ENDIF 02103 CALL dbcsr_insert_blocks (reduced, matrix, error=error) 02104 ! Now do the iterative folding in ceil(log_2(numnodes)) steps. 02105 nsteps = ceil_log2 (numnodes) 02106 ! Write folding info. 02107 !write(*,*)routineN//" ceil_log2", numnodes, nsteps,& 02108 ! 32, ceil_log2(32), 7, ceil_log2(7), 2, ceil_log2(2) 02109 !write(*,*)routineN//" Reduction target:", reduction_target 02110 most_spread = ISHFT(1, nsteps) 02111 my_norm_node = MODULO (mynode-reduction_target, numnodes) 02112 step = 1 02113 NULLIFY (mat_comm) 02114 log_steps: DO WHILE (step .LE. nsteps) 02115 ! 02116 half_spread = ISHFT (1, step-1) 02117 spread = 2 * half_spread 02118 ! Write step info 02119 !WRITE(*,*)routineN//" Step", step, half_spread, spread 02120 ! Check whether we participate in a send/recv according to the 02121 ! rules and according to the number of nodes. 02122 !--- 02123 i_recv = MODULO (my_norm_node, spread) .EQ. 0 02124 src_p = MODULO (mynode+half_spread, numnodes) 02125 i_recv = i_recv .AND. (my_norm_node+half_spread) .LT. numnodes 02126 i_send = MODULO (my_norm_node, spread) .EQ. half_spread 02127 dst_p = MODULO (mynode-half_spread, numnodes) 02128 i_send = i_send .AND. MODULO(my_norm_node,spread)-half_spread .GE. 0 02129 ! Write comm info 02130 !WRITE(*,*)routineN//" i_recv,send", i_recv, src_p, i_send, dst_p 02131 CALL dbcsr_assert ("NOT", i_recv .AND. i_send, dbcsr_fatal_level,& 02132 dbcsr_internal_error, routineN, & 02133 "Can not send and receive in the same step.", __LINE__, error=error) 02134 ! 02135 if_recv: IF (i_recv) THEN 02136 mat_comm => mat_recv 02137 CALL mp_recv (get_sizes, src_p, step, mp_group) 02138 CALL dbcsr_data_ensure_size (mat_comm%m%data_area, get_sizes(i_d),& 02139 nocopy=.TRUE., zero_pad=.FALSE., error=error) 02140 CALL ensure_array_size (mat_comm%m%index, lb=1, ub=get_sizes(i_i),& 02141 nocopy=.TRUE.,& 02142 memory_type=dbcsr_get_index_memory_type(mat_comm),& 02143 zero_pad=.TRUE., error=error) 02144 CALL mp_recv (mat_comm%m%index(1:get_sizes(i_i)), & 02145 src_p, step, mp_group) 02146 IF (get_sizes(i_d) .GT. 0) THEN 02147 CALL dbcsr_data_set_pointer (area=data_ptr,& 02148 rsize=get_sizes(i_d), csize=1,& 02149 pointee=mat_comm%m%data_area) 02150 CALL dbcsr_recv_any (data_ptr, src_p, tag=step, comm=mp_group,& 02151 error=error) 02152 ENDIF 02153 CALL dbcsr_repoint_index (mat_comm%m) 02154 mat_comm%m%valid = .TRUE. 02155 ! Now add the new data into the receiving matrix. 02156 CALL dbcsr_add_reserved (reduced, mat_comm, error=error) 02157 ENDIF if_recv 02158 if_send: IF (i_send) THEN 02159 mat_comm => reduced 02160 my_sizes(i_i) = SIZE (mat_comm%m%index) 02161 my_sizes(i_d) = dbcsr_get_data_size (mat_comm) 02162 CALL mp_send (my_sizes, dst_p, step, mp_group) 02163 CALL mp_send (mat_comm%m%index(1:my_sizes(i_i)), dst_p,& 02164 step, mp_group) 02165 IF (my_sizes(i_d) .GT. 0) THEN 02166 CALL dbcsr_data_set_pointer (area=data_ptr,& 02167 rsize=my_sizes(i_d), csize=1,& 02168 pointee=mat_comm%m%data_area) 02169 CALL dbcsr_send_any (data_ptr, dst_p, tag=step, comm=mp_group,& 02170 error=error) 02171 ENDIF 02172 ENDIF if_send 02173 step = step + 1 02174 ENDDO log_steps 02175 CALL dbcsr_data_clear_pointer (data_ptr) 02176 CALL dbcsr_data_release (data_ptr) 02177 CALL dbcsr_release (mat_recv) 02178 IF (debug_mod .AND. (mynode .EQ. reduction_target)) THEN 02179 cs2 = dbcsr_checksum (reduced, local=.TRUE., error=error) 02180 IF (ABS(cs1) .GT. comp_eps) THEN 02181 CALL dbcsr_assert (ABS((cs2 - cs1) / cs1) .LT. comp_eps,& 02182 dbcsr_warning_level, dbcsr_internal_error, routineN,& 02183 "Reduced checksum differs", __LINE__, error=error) 02184 ENDIF 02185 ENDIF 02186 CALL dbcsr_error_stop (error_handler, error) 02187 END SUBROUTINE dbcsr_reduce 02188 02189 02190 ! ***************************************************************************** 02196 SUBROUTINE dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c, error) 02197 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a, matrix_b 02198 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_c 02199 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02200 02201 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_hadamard_product', 02202 routineP = moduleN//':'//routineN 02203 02204 INTEGER :: blk, col, col_size, 02205 data_type, error_handler, 02206 nze, row, row_size 02207 LOGICAL :: found, tr_a, tr_b 02208 TYPE(dbcsr_data_obj) :: a_data, b_data, c_data 02209 TYPE(dbcsr_iterator) :: iter 02210 02211 ! --------------------------------------------------------------------------- 02212 02213 CALL dbcsr_error_set(routineN, error_handler, error) 02214 CALL dbcsr_assert (dbcsr_get_data_type(matrix_a).EQ.dbcsr_get_data_type(matrix_b).AND.& 02215 dbcsr_get_data_type(matrix_a).EQ.dbcsr_get_data_type(matrix_c), & 02216 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, & 02217 "data types not consistent, need to fix that",__LINE__,error) 02218 02219 CALL dbcsr_assert (dbcsr_nblkrows_total(matrix_a).EQ.dbcsr_nblkrows_total(matrix_b).AND.& 02220 dbcsr_nblkrows_total(matrix_c).EQ.dbcsr_nblkrows_total(matrix_a), & 02221 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, & 02222 "matrices not consistent",__LINE__,error) 02223 02224 data_type = dbcsr_get_data_type(matrix_a) 02225 CALL dbcsr_data_init (c_data) 02226 CALL dbcsr_data_new (c_data, data_type,& 02227 data_size=dbcsr_max_row_size(matrix_a)*dbcsr_max_col_size(matrix_a)) 02228 CALL dbcsr_set(matrix_c, dbcsr_scalar_zero(data_type), error=error) 02229 CALL dbcsr_data_init (a_data) 02230 CALL dbcsr_data_new (a_data, data_type) 02231 CALL dbcsr_data_init (b_data) 02232 CALL dbcsr_data_new (b_data, data_type) 02233 CALL dbcsr_iterator_start(iter, matrix_a) 02234 DO WHILE (dbcsr_iterator_blocks_left(iter)) 02235 SELECT CASE (dbcsr_get_data_type(matrix_a)) 02236 !CASE (dbcsr_type_real_4) 02237 CASE (dbcsr_type_real_8) 02238 CALL dbcsr_iterator_next_block(iter, row, col, a_data, tr_a, blk, & 02239 row_size=row_size, col_size=col_size) 02240 nze = row_size * col_size 02241 CALL dbcsr_get_block_p(matrix_b, row, col, b_data, tr_b, found) 02242 CALL dbcsr_assert (tr_a.EQV.tr_b, dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, & 02243 "tr not consistent, need to fix that",__LINE__,error) 02244 IF(found) THEN 02245 SELECT CASE (data_type) 02246 CASE (dbcsr_type_real_4) 02247 c_data%d%r_sp(1:nze) = a_data%d%r_sp(1:nze) * b_data%d%r_sp(1:nze) 02248 CASE (dbcsr_type_real_8) 02249 c_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze) * b_data%d%r_dp(1:nze) 02250 CASE (dbcsr_type_complex_4) 02251 c_data%d%c_sp(1:nze) = a_data%d%c_sp(1:nze) * b_data%d%c_sp(1:nze) 02252 CASE (dbcsr_type_complex_8) 02253 c_data%d%c_dp(1:nze) = a_data%d%c_dp(1:nze) * b_data%d%c_dp(1:nze) 02254 END SELECT 02255 CALL dbcsr_put_block(matrix_c, row, col, c_data, tr_a, & 02256 summation=.FALSE.) 02257 ENDIF 02258 !CASE (dbcsr_type_complex_4) 02259 !CASE (dbcsr_type_complex_8) 02260 CASE DEFAULT 02261 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 02262 routineN,"Only real double precision",__LINE__,error) 02263 END SELECT 02264 ENDDO 02265 CALL dbcsr_iterator_stop(iter) 02266 CALL dbcsr_finalize (matrix_c, error=error) 02267 CALL dbcsr_data_clear_pointer (a_data) 02268 CALL dbcsr_data_clear_pointer (b_data) 02269 CALL dbcsr_data_release (c_data) 02270 CALL dbcsr_data_release (a_data) 02271 CALL dbcsr_data_release (b_data) 02272 CALL dbcsr_error_stop(error_handler, error) 02273 END SUBROUTINE dbcsr_hadamard_product 02274 02275 ! ***************************************************************************** 02283 SUBROUTINE dbcsr_replace_blocks(matrix_a, matrix_b, error) 02284 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 02285 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 02286 TYPE(dbcsr_error_type), INTENT(inout) :: error 02287 02288 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_replace_blocks', 02289 routineP = moduleN//':'//routineN 02290 02291 COMPLEX(KIND=dp), DIMENSION(:, :), 02292 POINTER :: data_c 02293 COMPLEX(KIND=sp), DIMENSION(:, :), 02294 POINTER :: data_z 02295 INTEGER :: blk, col, error_handler, row 02296 LOGICAL :: tr 02297 REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_d 02298 REAL(KIND=sp), DIMENSION(:, :), POINTER :: data_r 02299 TYPE(dbcsr_iterator) :: iter 02300 02301 ! --------------------------------------------------------------------------- 02302 02303 CALL dbcsr_error_set(routineN, error_handler, error) 02304 02305 CALL dbcsr_assert (dbcsr_get_data_type(matrix_a).EQ.dbcsr_get_data_type(matrix_b), & 02306 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, "data types not consistent",__LINE__,error) 02307 ! 02308 ! let's go 02309 CALL dbcsr_assert (matrix_a%m%nblkrows_total.EQ.matrix_b%m%nblkrows_total, & 02310 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, "matrices not consistent",__LINE__,error) 02311 02312 CALL dbcsr_iterator_start(iter, matrix_b) 02313 02314 DO WHILE (dbcsr_iterator_blocks_left(iter)) 02315 02316 SELECT CASE (dbcsr_get_data_type(matrix_a)) 02317 CASE (dbcsr_type_real_4) 02318 CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk) 02319 CALL dbcsr_put_block(matrix_a, row, col, data_r) 02320 CASE (dbcsr_type_real_8) 02321 CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk) 02322 CALL dbcsr_put_block(matrix_a, row, col, data_d, tr,& 02323 summation=.FALSE.) 02324 CASE (dbcsr_type_complex_4) 02325 CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk) 02326 CALL dbcsr_put_block(matrix_a, row, col, data_c) 02327 CASE (dbcsr_type_complex_8) 02328 CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk) 02329 CALL dbcsr_put_block(matrix_a, row, col, data_z) 02330 CASE DEFAULT 02331 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 02332 routineN,"Unkown precision",__LINE__,error) 02333 END SELECT 02334 02335 ENDDO 02336 02337 CALL dbcsr_iterator_stop(iter) 02338 CALL dbcsr_finalize (matrix_a, error=error) 02339 ! 02340 CALL dbcsr_error_stop(error_handler, error) 02341 END SUBROUTINE dbcsr_replace_blocks 02342 02343 ! ***************************************************************************** 02348 SUBROUTINE dbcsr_add_on_diag_anytype(matrix, alpha_scalar, first_row, last_row, error) 02349 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 02350 TYPE(dbcsr_scalar_type), INTENT(IN) :: alpha_scalar 02351 INTEGER, INTENT(in), OPTIONAL :: first_row, last_row 02352 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02353 02354 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag_anytype', 02355 routineP = moduleN//':'//routineN 02356 02357 INTEGER :: error_handler, hold, imax, imin, my_first_row, my_last_row, 02358 mynode, offset_beg, offset_end, row, row_size, stored_row 02359 INTEGER, DIMENSION(:), POINTER :: row_blk_offsets 02360 LOGICAL :: found, tr 02361 TYPE(dbcsr_data_obj) :: buff, data_a, small_buff 02362 02363 ! --------------------------------------------------------------------------- 02364 02365 CALL dbcsr_error_set(routineN, error_handler, error) 02366 CALL dbcsr_assert (dbcsr_nblkrows_total(matrix).EQ.dbcsr_nblkcols_total(matrix).AND.& 02367 dbcsr_nfullrows_total(matrix).EQ.dbcsr_nfullrows_total(matrix), & 02368 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, "matrices not consistent",__LINE__,error) 02369 02370 my_first_row = 1 02371 my_last_row = dbcsr_nfullrows_total(matrix) 02372 IF(PRESENT(first_row)) my_first_row = first_row 02373 IF(PRESENT(last_row)) my_last_row = last_row 02374 02375 mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) 02376 row_blk_offsets => array_data(dbcsr_col_block_offsets(matrix)) 02377 02378 CALL dbcsr_work_create(matrix, work_mutable=.TRUE., error=error) 02379 CALL dbcsr_data_init (buff) 02380 CALL dbcsr_data_init (data_a) 02381 CALL dbcsr_data_init (small_buff) 02382 CALL dbcsr_data_new (data_a,& 02383 dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix))) 02384 CALL dbcsr_data_new (buff,& 02385 dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix)),& 02386 dbcsr_max_row_size(matrix), dbcsr_max_col_size(matrix)) 02387 CALL dbcsr_data_new (small_buff,& 02388 dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix))) 02389 02390 DO row = 1,dbcsr_nblkrows_total(matrix) 02391 tr = .FALSE. 02392 stored_row = row 02393 CALL dbcsr_get_stored_coordinates (matrix, stored_row, stored_row, tr,& 02394 hold) 02395 IF(hold.EQ.mynode) THEN 02396 CALL dbcsr_get_block_p(matrix, stored_row, stored_row, data_a, tr,& 02397 found, row_size=row_size) 02398 offset_beg = row_blk_offsets(row) 02399 offset_end = row_blk_offsets(row+1)-1 02400 IF(my_first_row.GT.offset_end.OR.my_last_row.LT.offset_beg) CYCLE 02401 imin = 1 02402 IF(my_first_row.gt.offset_beg) THEN 02403 imin = my_first_row - offset_beg + 1 02404 ENDIF 02405 imax = row_size 02406 IF(my_last_row.lt.offset_end) THEN 02407 imax = my_last_row - offset_end + row_size 02408 ENDIF 02409 IF(found) THEN 02410 CALL block_add_on_diag(data_a, alpha_scalar, row_size, & 02411 imin=imin, imax=imax, error=error) 02412 ELSE 02413 CALL dbcsr_data_set_pointer (small_buff, row_size, row_size,& 02414 buff) 02415 CALL dbcsr_data_clear (small_buff) 02416 CALL block_add_on_diag (small_buff, alpha_scalar, row_size,& 02417 imin=imin, imax=imax, error=error) 02418 CALL dbcsr_put_block (matrix, stored_row, stored_row, small_buff) 02419 ENDIF 02420 ENDIF 02421 ENDDO 02422 02423 CALL dbcsr_data_clear_pointer (data_a) 02424 CALL dbcsr_data_clear_pointer (small_buff) 02425 CALL dbcsr_data_release (small_buff) 02426 CALL dbcsr_data_release (buff) 02427 CALL dbcsr_data_release (data_a) 02428 CALL dbcsr_finalize(matrix, error=error) 02429 CALL dbcsr_error_stop(error_handler, error) 02430 END SUBROUTINE dbcsr_add_on_diag_anytype 02431 02432 SUBROUTINE dbcsr_init_random(matrix, error) 02433 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 02434 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02435 02436 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_init_random', 02437 routineP = moduleN//':'//routineN 02438 02439 INTEGER :: col, col_size, error_handler, 02440 hold, iseed(4), mynode, row, 02441 row_size, stored_col, 02442 stored_row 02443 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size 02444 LOGICAL :: found, tr 02445 REAL(real_8), ALLOCATABLE, DIMENSION(:) :: rnd 02446 REAL(real_8), DIMENSION(:, :), POINTER :: buff, data_d 02447 02448 ! --------------------------------------------------------------------------- 02449 02450 CALL dbcsr_error_set(routineN, error_handler, error) 02451 02452 row_blk_size => array_data (matrix%m%row_blk_size) 02453 col_blk_size => array_data (matrix%m%col_blk_size) 02454 mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) 02455 CALL dbcsr_work_create(matrix, work_mutable=.TRUE., error=error) 02456 02457 iseed(1)=4;iseed(2)=3;iseed(3)=2;iseed(4)=1! set the seed for dlarnv 02458 ALLOCATE(rnd(MAXVAL(row_blk_size)*MAXVAL(col_blk_size))) 02459 DO row = 1,dbcsr_nblkrows_total(matrix) 02460 DO col = 1,dbcsr_nblkcols_total(matrix) 02461 row_size = row_blk_size(row) 02462 col_size = col_blk_size(col) 02463 CALL dlarnv(1,iseed,row_size*col_size,rnd(1)) 02464 tr = .FALSE. 02465 stored_row = row 02466 stored_col = col 02467 CALL dbcsr_get_stored_coordinates(matrix, stored_row, stored_col, tr, hold) 02468 IF(hold.EQ.mynode) THEN 02469 CALL dbcsr_get_block_p(matrix, stored_row, stored_col, data_d, tr, found) 02470 IF(found) THEN 02471 CALL dcopy(row_size*col_size,rnd,1,data_d,1) 02472 ELSE 02473 ALLOCATE(buff(row_size,col_size)) 02474 CALL dcopy(row_size*col_size,rnd,1,buff,1) 02475 CALL dbcsr_put_block (matrix, stored_row, stored_col, buff) 02476 DEALLOCATE(buff) 02477 ENDIF 02478 ENDIF 02479 ENDDO 02480 ENDDO 02481 DEALLOCATE(rnd) 02482 02483 CALL dbcsr_finalize(matrix, error=error) 02484 CALL dbcsr_error_stop(error_handler, error) 02485 02486 END SUBROUTINE dbcsr_init_random 02487 02488 ! ***************************************************************************** 02494 SUBROUTINE dbcsr_get_block_diag(matrix, diag, error) 02495 02496 TYPE(dbcsr_obj), INTENT(IN) :: matrix 02497 TYPE(dbcsr_obj), INTENT(INOUT) :: diag 02498 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02499 02500 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_block_diag', 02501 routineP = moduleN//':'//routineN 02502 02503 INTEGER :: blk, col, error_handler, row 02504 LOGICAL :: tr 02505 TYPE(dbcsr_data_obj) :: data_a 02506 TYPE(dbcsr_iterator) :: iter 02507 02508 ! --------------------------------------------------------------------------- 02509 02510 CALL dbcsr_error_set(routineN, error_handler, error) 02511 CALL dbcsr_create(diag, name='diag of '//TRIM(matrix%m%name), & 02512 template=matrix, error=error) 02513 02514 CALL dbcsr_data_init (data_a) 02515 CALL dbcsr_data_new (data_a, dbcsr_get_data_type(matrix)) 02516 CALL dbcsr_iterator_start(iter, matrix) 02517 DO WHILE (dbcsr_iterator_blocks_left(iter)) 02518 CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr, blk) 02519 IF(row.EQ.col) CALL dbcsr_put_block(diag, row, col, data_a, tr) 02520 ENDDO 02521 CALL dbcsr_iterator_stop(iter) 02522 CALL dbcsr_data_clear_pointer (data_a) 02523 CALL dbcsr_data_release (data_a) 02524 CALL dbcsr_finalize(diag, error=error) 02525 CALL dbcsr_error_stop(error_handler, error) 02526 END SUBROUTINE dbcsr_get_block_diag 02527 02528 02529 ! ***************************************************************************** 02535 SUBROUTINE dbcsr_get_diag_anytype(matrix, diag, error) 02536 02537 TYPE(dbcsr_obj), INTENT(IN) :: matrix 02538 TYPE(dbcsr_data_obj), INTENT(INOUT) :: diag 02539 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02540 02541 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag_anytype', 02542 routineP = moduleN//':'//routineN 02543 02544 INTEGER :: blk, col, data_type, 02545 error_handler, row, 02546 row_offset, row_size 02547 LOGICAL :: tr 02548 TYPE(dbcsr_data_obj) :: data_a, diag_a 02549 TYPE(dbcsr_iterator) :: iter 02550 02551 ! --------------------------------------------------------------------------- 02552 02553 CALL dbcsr_error_set(routineN, error_handler, error) 02554 02555 CALL dbcsr_assert (dbcsr_nfullrows_total(matrix), "LE", dbcsr_data_get_size(diag), & 02556 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02557 "Diagonal is too small",__LINE__,error) 02558 02559 CALL dbcsr_data_clear (diag) 02560 data_type = dbcsr_get_data_type (matrix) 02561 02562 CALL dbcsr_data_init (data_a) 02563 CALL dbcsr_data_new (data_a,& 02564 dbcsr_type_1d_to_2d (data_type)) 02565 CALL dbcsr_data_init (diag_a) 02566 CALL dbcsr_data_new (diag_a, data_type) 02567 02568 CALL dbcsr_iterator_start(iter, matrix) 02569 DO WHILE (dbcsr_iterator_blocks_left(iter)) 02570 CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr, blk, & 02571 row_offset=row_offset, row_size=row_size) 02572 IF(row.EQ.col) THEN 02573 diag_a = pointer_view (diag_a, diag, offset=row_offset, len=row_size) 02574 CALL get_block2d_diagonal (data_a, diag_a, row_size, error=error) 02575 ENDIF 02576 ENDDO 02577 CALL dbcsr_iterator_stop(iter) 02578 02579 CALL dbcsr_data_clear_pointer (diag_a) 02580 CALL dbcsr_data_release (diag_a) 02581 CALL dbcsr_data_clear_pointer (data_a) 02582 CALL dbcsr_data_release (data_a) 02583 CALL dbcsr_error_stop(error_handler, error) 02584 END SUBROUTINE dbcsr_get_diag_anytype 02585 02586 ! ***************************************************************************** 02592 SUBROUTINE dbcsr_set_diag_anytype(matrix, diag, error) 02593 02594 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 02595 TYPE(dbcsr_data_obj), INTENT(IN) :: diag 02596 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02597 02598 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag_anytype', 02599 routineP = moduleN//':'//routineN 02600 02601 INTEGER :: blk, col, dt, error_handler, 02602 row, row_offset, row_size 02603 LOGICAL :: tr 02604 TYPE(dbcsr_data_obj) :: data_a, diag_a 02605 TYPE(dbcsr_iterator) :: iter 02606 02607 ! --------------------------------------------------------------------------- 02608 02609 CALL dbcsr_error_set(routineN, error_handler, error) 02610 02611 CALL dbcsr_assert (dbcsr_nfullrows_total(matrix), "LE", dbcsr_data_get_size(diag), & 02612 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02613 "Diagonal too small",__LINE__,error) 02614 02615 dt = dbcsr_get_data_type (matrix) 02616 CALL dbcsr_data_init (data_a) 02617 CALL dbcsr_data_new (data_a, dbcsr_type_1d_to_2d (dt)) 02618 CALL dbcsr_data_init (diag_a) 02619 CALL dbcsr_data_new (diag_a, dt) 02620 02621 CALL dbcsr_iterator_start(iter, matrix) 02622 DO WHILE (dbcsr_iterator_blocks_left(iter)) 02623 CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr, blk, & 02624 row_offset=row_offset, row_size=row_size) 02625 IF(row.EQ.col) THEN 02626 diag_a = pointer_view (diag_a, diag, offset=row_offset, len=row_size) 02627 CALL set_block2d_diagonal (data_a, diag_a, row_size, error=error) 02628 ENDIF 02629 ENDDO 02630 CALL dbcsr_iterator_stop(iter) 02631 02632 CALL dbcsr_data_clear_pointer (diag_a) 02633 CALL dbcsr_data_release (diag_a) 02634 CALL dbcsr_data_clear_pointer (data_a) 02635 CALL dbcsr_data_release (data_a) 02636 CALL dbcsr_error_stop(error_handler, error) 02637 END SUBROUTINE dbcsr_set_diag_anytype 02638 02639 02640 ! ***************************************************************************** 02644 LOGICAL FUNCTION symmetry_consistent(matrix_type,data_type,error) 02645 CHARACTER, INTENT(IN) :: matrix_type 02646 INTEGER, INTENT(IN) :: data_type 02647 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02648 02649 CHARACTER(len=*), PARAMETER :: routineN = 'symmetry_consistent', 02650 routineP = moduleN//':'//routineN 02651 02652 symmetry_consistent = .FALSE. 02653 02654 SELECT CASE (data_type) 02655 CASE(dbcsr_type_real_4,dbcsr_type_real_8) 02656 SELECT CASE (matrix_type) 02657 CASE (dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_type_antisymmetric) 02658 symmetry_consistent = .TRUE. 02659 END SELECT 02660 CASE(dbcsr_type_complex_4,dbcsr_type_complex_8) 02661 SELECT CASE (matrix_type) 02662 CASE (dbcsr_type_no_symmetry, dbcsr_type_hermitian, dbcsr_type_antihermitian) 02663 symmetry_consistent = .TRUE. 02664 END SELECT 02665 CASE DEFAULT 02666 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error,& 02667 routineN, "Invalid data type.",__LINE__,error) 02668 END SELECT 02669 02670 END FUNCTION symmetry_consistent 02671 02672 ! ***************************************************************************** 02676 LOGICAL FUNCTION symmetry_compatible(matrix_type_a,matrix_type_b,error) 02677 CHARACTER, INTENT(IN) :: matrix_type_a, matrix_type_b 02678 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02679 02680 CHARACTER(len=*), PARAMETER :: routineN = 'symmetry_compatible', 02681 routineP = moduleN//':'//routineN 02682 02683 symmetry_compatible = .FALSE. 02684 02685 SELECT CASE (matrix_type_a) 02686 CASE (dbcsr_type_no_symmetry) 02687 SELECT CASE(matrix_type_b) 02688 CASE(dbcsr_type_no_symmetry) 02689 symmetry_compatible = .TRUE. 02690 END SELECT 02691 CASE(dbcsr_type_symmetric, dbcsr_type_hermitian) 02692 SELECT CASE(matrix_type_b) 02693 CASE(dbcsr_type_symmetric, dbcsr_type_hermitian) 02694 symmetry_compatible = .TRUE. 02695 END SELECT 02696 CASE(dbcsr_type_antisymmetric, dbcsr_type_antihermitian) 02697 SELECT CASE(matrix_type_b) 02698 CASE(dbcsr_type_antisymmetric, dbcsr_type_antihermitian) 02699 symmetry_compatible = .TRUE. 02700 END SELECT 02701 CASE DEFAULT 02702 CALL dbcsr_assert(.FALSE., dbcsr_failure_level,& 02703 dbcsr_wrong_args_error, routineP, "Invalid matrix type.",__LINE__,error) 02704 END SELECT 02705 02706 END FUNCTION symmetry_compatible 02707 02708 02709 ! ***************************************************************************** 02724 SUBROUTINE dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity,& 02725 shallow_data, keep_imaginary, matrix_type, error) 02726 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 02727 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 02728 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name 02729 LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity, shallow_data, 02730 keep_imaginary 02731 CHARACTER, INTENT(IN), OPTIONAL :: matrix_type 02732 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02733 02734 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy', 02735 routineP = moduleN//':'//routineN 02736 02737 CHARACTER :: new_matrix_type, repl_type 02738 INTEGER :: error_handler, new_type 02739 LOGICAL :: keep_sparse, shallow 02740 02741 ! --------------------------------------------------------------------------- 02742 02743 CALL dbcsr_error_set(routineN, error_handler, error) 02744 CALL dbcsr_assert (symmetry_consistent(dbcsr_get_matrix_type(matrix_a),& 02745 dbcsr_get_data_type(matrix_a), error),& 02746 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 02747 "Source matrix symmetry not consistent with its data type.",__LINE__,error) 02748 shallow = .FALSE. ; IF (PRESENT (shallow_data)) shallow = shallow_data 02749 keep_sparse = .FALSE. 02750 IF (PRESENT (keep_sparsity)) keep_sparse = keep_sparsity 02751 CALL dbcsr_access_flush (matrix_a, error=error) 02752 !CALL dbcsr_assert (dbcsr_is_initialized (matrix_b),& 02753 ! dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02754 ! "Target matrix must be initialized", error=error) 02755 CALL dbcsr_assert (.not.keep_sparse.or.dbcsr_valid_index(matrix_b),& 02756 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02757 "Target matrix must be valid to keep its sparsity",__LINE__,error) 02758 CALL dbcsr_assert (.not.keep_sparse.or..not.shallow, dbcsr_warning_level,& 02759 dbcsr_wrong_args_error, routineN,& 02760 "Shallow copy not compatibly with sparsity retainment",__LINE__,error) 02761 IF (keep_sparse) THEN 02762 IF (PRESENT (name)) matrix_b%m%name = name 02763 CALL dbcsr_copy_into_existing(matrix_b, matrix_a, error) 02764 ELSE 02765 IF (dbcsr_is_initialized (matrix_b)) THEN 02766 new_type = dbcsr_get_data_type (matrix_b) 02767 repl_type = dbcsr_get_replication_type(matrix_b) 02768 ELSE 02769 new_type = dbcsr_get_data_type (matrix_a) 02770 repl_type = dbcsr_get_replication_type(matrix_a) 02771 ENDIF 02772 new_matrix_type = dbcsr_get_matrix_type (matrix_a) 02773 IF (PRESENT (matrix_type)) THEN 02774 CALL dbcsr_assert (symmetry_compatible(dbcsr_get_matrix_type(matrix_a),& 02775 matrix_type, error),& 02776 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02777 "Specified target matrix symmetry "//matrix_type//& 02778 " not compatible with source matrix type "//dbcsr_get_matrix_type(matrix_a),& 02779 __LINE__,error) 02780 new_matrix_type = matrix_type 02781 END IF 02782 CALL dbcsr_assert (symmetry_consistent(new_matrix_type,new_type,error),& 02783 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02784 "Target matrix symmetry "//new_matrix_type//" not consistent with its data type.",& 02785 __LINE__,error) 02786 IF(PRESENT(name)) THEN 02787 CALL dbcsr_create(matrix_b, name=TRIM(name), & 02788 template = matrix_a,& 02789 matrix_type = new_matrix_type,& 02790 data_type = new_type,& 02791 error=error) 02792 ELSE 02793 CALL dbcsr_create(matrix_b,& 02794 name='copy of '//TRIM(dbcsr_name(matrix_a)),& 02795 data_type = new_type,& 02796 matrix_type = new_matrix_type,& 02797 template = matrix_a, error=error) 02798 ENDIF 02799 CALL ensure_array_size(matrix_b%m%index, ub=SIZE(matrix_a%m%index),& 02800 memory_type = dbcsr_get_index_memory_type(matrix_b), error=error) 02801 ! 02802 ! copy index and data 02803 matrix_b%m%index(1:SIZE(matrix_a%m%index)) = matrix_a%m%index(:) 02804 IF (.NOT. shallow) THEN 02805 CALL dbcsr_assert (matrix_a%m%nze, "LE", dbcsr_get_data_size(matrix_a),& 02806 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 02807 "Source matrix sizes not consistent!",__LINE__,error) 02808 CALL dbcsr_data_ensure_size (matrix_b%m%data_area,& 02809 dbcsr_get_data_size(matrix_a), error=error) 02810 IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_get_data_type(matrix_b))& 02811 THEN 02812 CALL dbcsr_data_copyall (matrix_b%m%data_area,& 02813 matrix_a%m%data_area, error=error) 02814 ELSE 02815 CALL dbcsr_data_convert (matrix_b%m%data_area,& 02816 matrix_a%m%data_area, drop_real=keep_imaginary) 02817 ENDIF 02818 ELSE 02819 CALL dbcsr_assert (dbcsr_get_data_type(matrix_a) & 02820 .EQ. dbcsr_get_data_type(matrix_b), dbcsr_fatal_level,& 02821 dbcsr_wrong_args_error, routineN,& 02822 "Shallow copy only possible when retaining data type.", __LINE__, error) 02823 CALL dbcsr_switch_data_area (matrix_b, matrix_a%m%data_area,& 02824 error=error) 02825 ENDIF 02826 ! 02827 ! the row_p, col_i and blk_p ... 02828 CALL dbcsr_repoint_index(matrix_b%m) 02829 matrix_b%m%nze = matrix_a%m%nze 02830 matrix_b%m%nblks = matrix_b%m%nblks 02831 matrix_b%m%valid = .TRUE. 02832 02833 matrix_b%m%sparsity_id = matrix_a%m%sparsity_id 02834 ENDIF 02835 CALL dbcsr_error_stop(error_handler, error) 02836 END SUBROUTINE dbcsr_copy 02837 02838 ! ***************************************************************************** 02844 SUBROUTINE dbcsr_copy_into_existing(matrix_b, matrix_a, error) 02845 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 02846 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 02847 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02848 02849 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_into_existing', 02850 routineP = moduleN//':'//routineN 02851 02852 INTEGER :: col_size, data_type, dst_col, dst_row, error_handler, rel, 02853 row_size, src_col, src_cs, src_row, src_rs 02854 LOGICAL :: dst_tr, making_symmetric, 02855 neg_imag, neg_real, src_tr 02856 TYPE(dbcsr_data_obj) :: dst_data, src_data 02857 TYPE(dbcsr_iterator) :: dst_iter, src_iter 02858 02859 ! --------------------------------------------------------------------------- 02860 02861 CALL dbcsr_error_set(routineN, error_handler, error) 02862 CALL dbcsr_assert (dbcsr_get_data_type(matrix_b)& 02863 .EQ. dbcsr_get_data_type(matrix_a), dbcsr_fatal_level,& 02864 dbcsr_wrong_args_error, routineN, "Matrices have different data types.",__LINE__,error) 02865 data_type = dbcsr_get_data_type (matrix_b) 02866 neg_real = matrix_b%m%negate_real 02867 neg_imag = matrix_b%m%negate_imaginary 02868 making_symmetric = dbcsr_has_symmetry (matrix_b)& 02869 .AND. .NOT. dbcsr_has_symmetry (matrix_a) 02870 IF (making_symmetric) THEN 02871 CALL dbcsr_copy_into_existing_sym(matrix_b, matrix_a, error) 02872 CALL dbcsr_error_stop(error_handler, error) 02873 RETURN 02874 ENDIF 02875 CALL dbcsr_data_init (src_data) 02876 CALL dbcsr_data_init (dst_data) 02877 CALL dbcsr_data_new (src_data, data_type) 02878 CALL dbcsr_data_new (dst_data, data_type) 02879 CALL dbcsr_iterator_start (src_iter, matrix_a) 02880 CALL dbcsr_iterator_start (dst_iter, matrix_b) 02881 ! Iterate through the blocks of the source and destination 02882 ! matrix. There are three possibilites: 1. copy the data for 02883 ! blocks present in both; 2 skip source blocks not present in the 02884 ! target; 3 zero blocks not present in the source. 02885 IF (dbcsr_iterator_blocks_left (src_iter)) THEN 02886 CALL dbcsr_iterator_next_block (src_iter, src_row, src_col, src_data,& 02887 src_tr) 02888 ELSE 02889 src_row = 0 ; src_col = 0 02890 ENDIF 02891 DO WHILE (dbcsr_iterator_blocks_left (dst_iter)) 02892 CALL dbcsr_iterator_next_block (dst_iter, dst_row, dst_col, dst_data,& 02893 dst_tr, row_size=row_size, col_size=col_size) 02894 ! Now find the source position that is greater or equal to the 02895 ! target one. I.e, skip blocks that the target doesn't have. 02896 rel = pos_relation (dst_row, dst_col, src_row, src_col) 02897 DO WHILE (rel .EQ. 1 .AND. dbcsr_iterator_blocks_left (src_iter)) 02898 CALL dbcsr_iterator_next_block (src_iter, src_row, src_col,& 02899 src_data, src_tr, row_size=src_rs, col_size=src_cs) 02900 rel = pos_relation (dst_row, dst_col, src_row, src_col) 02901 ENDDO 02902 SELECT CASE (rel) 02903 CASE (-1, 1) 02904 ! Target lags source or ran out of source 02905 CALL dbcsr_data_clear (dst_data) 02906 CASE (0) 02907 ! Copy the data 02908 CALL dbcsr_assert (dbcsr_data_get_size (src_data)& 02909 .EQ. dbcsr_data_get_size (dst_data), dbcsr_fatal_level,& 02910 dbcsr_internal_error, routineN, "Block sizes not equal!",__LINE__,error) 02911 IF (src_tr .EQV. dst_tr) THEN 02912 CALL dbcsr_data_copyall (dst_data, src_data, error=error) 02913 ELSE 02914 CALL dbcsr_block_partial_copy(dst=dst_data, dst_tr = dst_tr,& 02915 dst_rs = row_size, dst_cs = col_size,& 02916 dst_r_lb = 1, dst_c_lb = 1,& 02917 src = src_data, src_tr = src_tr,& 02918 src_rs = src_rs, src_cs = src_cs,& 02919 src_r_lb = 1, src_c_lb = 1,& 02920 nrow = row_size, ncol = col_size) 02921 IF (neg_real) THEN 02922 CALL dbcsr_block_real_neg (dst_data, row_size, col_size, error=error) 02923 ENDIF 02924 IF (neg_imag) THEN 02925 CALL dbcsr_block_conjg (dst_data, row_size, col_size, error=error) 02926 ENDIF 02927 ENDIF 02928 CASE default 02929 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_internal_error,& 02930 routineN, "Trouble syncing iterators",__LINE__,error) 02931 END SELECT 02932 END DO 02933 CALL dbcsr_iterator_stop (src_iter) 02934 CALL dbcsr_iterator_stop (dst_iter) 02935 CALL dbcsr_data_clear_pointer (src_data) 02936 CALL dbcsr_data_clear_pointer (dst_data) 02937 CALL dbcsr_data_release (src_data) 02938 CALL dbcsr_data_release (dst_data) 02939 CALL dbcsr_error_stop(error_handler, error) 02940 02941 END SUBROUTINE dbcsr_copy_into_existing 02942 02943 ! ***************************************************************************** 02949 SUBROUTINE dbcsr_copy_into_existing_sym(matrix_b, matrix_a, error) 02950 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 02951 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 02952 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 02953 02954 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_into_existing_sym', 02955 routineP = moduleN//':'//routineN 02956 02957 INTEGER :: col_size, data_type, dst_col, dst_row, error_handler, 02958 row_size, src_col, src_cs, src_row, src_rs 02959 LOGICAL :: dst_tr, found, neg_imag, 02960 neg_real, src_tr 02961 TYPE(dbcsr_data_obj) :: dst_data, src_data 02962 TYPE(dbcsr_iterator) :: dst_iter 02963 02964 ! --------------------------------------------------------------------------- 02965 02966 CALL dbcsr_error_set(routineN, error_handler, error) 02967 CALL dbcsr_assert (dbcsr_get_data_type(matrix_b)& 02968 .EQ. dbcsr_get_data_type(matrix_a), dbcsr_fatal_level,& 02969 dbcsr_wrong_args_error, routineN, "Matrices have different data types.",__LINE__,error) 02970 data_type = dbcsr_get_data_type (matrix_b) 02971 CALL dbcsr_assert (dbcsr_has_symmetry (matrix_b)& 02972 .AND. .NOT. dbcsr_has_symmetry (matrix_a),& 02973 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 02974 "Must copy from non-symmetric to symmetric matrix.",& 02975 __LINE__, error=error) 02976 neg_real = matrix_b%m%negate_real 02977 neg_imag = matrix_b%m%negate_imaginary 02978 02979 CALL dbcsr_data_init (src_data) 02980 CALL dbcsr_data_init (dst_data) 02981 CALL dbcsr_data_new (src_data, data_type) 02982 CALL dbcsr_data_new (dst_data, data_type) 02983 CALL dbcsr_iterator_start (dst_iter, matrix_b) 02984 ! Iterate through the blocks of the destination matrix. For each 02985 ! one, try to find an appropriate source matrix block and copy it 02986 ! into the destination matrix. 02987 DO WHILE (dbcsr_iterator_blocks_left (dst_iter)) 02988 CALL dbcsr_iterator_next_block (dst_iter, dst_row, dst_col, dst_data,& 02989 dst_tr, row_size=row_size, col_size=col_size) 02990 src_row = dst_row 02991 src_col = dst_col 02992 IF (checker_tr (dst_row, dst_col))& 02993 CALL swap (src_row, src_col) 02994 CALL dbcsr_get_block_p (matrix_a, src_row, src_col, src_data, src_tr,& 02995 found=found, row_size=src_rs, col_size=src_cs) 02996 IF (.NOT. found) THEN 02997 CALL dbcsr_data_clear (dst_data) 02998 ELSE 02999 IF (checker_tr (dst_row, dst_col)) THEN 03000 src_tr = .NOT. src_tr 03001 CALL swap (src_rs, src_cs) 03002 ENDIF 03003 CALL dbcsr_block_partial_copy(dst=dst_data, dst_tr = dst_tr,& 03004 dst_rs = row_size, dst_cs = col_size,& 03005 dst_r_lb = 1, dst_c_lb = 1,& 03006 src = src_data, src_tr = src_tr,& 03007 src_rs = src_rs, src_cs = src_cs,& 03008 src_r_lb = 1, src_c_lb = 1,& 03009 nrow = row_size, ncol = col_size) 03010 IF (neg_real) THEN 03011 CALL dbcsr_block_real_neg (dst_data, row_size, col_size, error=error) 03012 ENDIF 03013 IF (neg_imag) THEN 03014 CALL dbcsr_block_conjg (dst_data, row_size, col_size, error=error) 03015 ENDIF 03016 ENDIF 03017 END DO 03018 CALL dbcsr_iterator_stop (dst_iter) 03019 CALL dbcsr_data_clear_pointer (src_data) 03020 CALL dbcsr_data_clear_pointer (dst_data) 03021 CALL dbcsr_data_release (src_data) 03022 CALL dbcsr_data_release (dst_data) 03023 CALL dbcsr_error_stop(error_handler, error) 03024 03025 END SUBROUTINE dbcsr_copy_into_existing_sym 03026 03027 03028 03029 ! ***************************************************************************** 03042 SUBROUTINE dbcsr_copy_columns(matrix_b, matrix_a,& 03043 ncol, source_start, target_start, error) 03044 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 03045 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 03046 INTEGER, INTENT(IN) :: ncol, source_start, 03047 target_start 03048 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03049 03050 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_columns', 03051 routineP = moduleN//':'//routineN 03052 03053 INTEGER :: col_size, data_type, dst_col, dst_row, error_handler, 03054 first_col, last_col, num_col, rel, row_size, src_col, src_row 03055 INTEGER, DIMENSION(:), POINTER :: dst_col_offsets, 03056 dst_col_sizes, 03057 src_col_offsets, src_col_sizes 03058 LOGICAL :: dst_tr, src_tr 03059 TYPE(dbcsr_data_obj) :: dst_data, src_data, tmp_buffer 03060 TYPE(dbcsr_iterator) :: dst_iter, src_iter 03061 03062 ! --------------------------------------------------------------------------- 03063 03064 CALL dbcsr_error_set(routineN, error_handler, error) 03065 CALL dbcsr_assert (dbcsr_get_data_type(matrix_b)& 03066 .EQ. dbcsr_get_data_type(matrix_a), dbcsr_fatal_level,& 03067 dbcsr_wrong_args_error, routineN, "Matrices have different data types.",__LINE__,error) 03068 data_type = dbcsr_get_data_type (matrix_b) 03069 03070 CALL dbcsr_assert (source_start .EQ. target_start, dbcsr_fatal_level,& 03071 dbcsr_unimplemented_error_nr, routineN, "Column shifting not supported",__LINE__,error) 03072 03073 src_col_offsets => array_data (matrix_a%m%col_blk_offset) 03074 dst_col_offsets => array_data (matrix_b%m%col_blk_offset) 03075 src_col_sizes => array_data (matrix_a%m%col_blk_size) 03076 dst_col_sizes => array_data (matrix_b%m%col_blk_size) 03077 03078 CALL dbcsr_data_init (src_data) 03079 CALL dbcsr_data_init (dst_data) 03080 CALL dbcsr_data_init (tmp_buffer) 03081 CALL dbcsr_data_new (src_data, data_type) 03082 CALL dbcsr_data_new (dst_data, data_type) 03083 CALL dbcsr_data_new (tmp_buffer, data_type,& 03084 matrix_a%m%max_rbs*matrix_a%m%max_cbs) 03085 CALL dbcsr_iterator_start (src_iter, matrix_a) 03086 CALL dbcsr_iterator_start (dst_iter, matrix_b) 03087 ! Iterate through the blocks of the source and destination 03088 ! matrix. There are three possibilites: 1 copy data for blocks 03089 ! present in both that are fully within the specified range; 2 03090 ! copy partial data for blacks partially within the specified 03091 ! range; 3 add (partial or full) data for blocks present in source 03092 ! but not in the target 03093 IF (dbcsr_iterator_blocks_left (src_iter)) THEN 03094 CALL dbcsr_iterator_next_block (src_iter, src_row, src_col, src_data,& 03095 src_tr) 03096 ELSE 03097 src_row = 0 ; src_col = 0 03098 ENDIF 03099 DO WHILE (dbcsr_iterator_blocks_left (dst_iter)) 03100 CALL dbcsr_iterator_next_block (dst_iter, dst_row, dst_col, dst_data,& 03101 dst_tr, row_size=row_size, col_size=col_size) 03102 ! Now find the source position that is greater or equal to the 03103 ! target one. I.e, skip blocks that the target doesn't have. 03104 rel = pos_relation (dst_row, dst_col, src_row, src_col) 03105 DO WHILE (rel .EQ. 1 .AND. dbcsr_iterator_blocks_left (src_iter)) 03106 CALL dbcsr_iterator_next_block (src_iter, src_row, src_col,& 03107 src_data, src_tr) 03108 rel = pos_relation (dst_row, dst_col, src_row, src_col) 03109 ENDDO 03110 SELECT CASE (rel) 03111 CASE (0) ! Equal block coordinates; block exists in both src and dst 03112 ! Check if the full block is being copied 03113 CALL dbcsr_assert (src_col_sizes(src_col) .EQ. dst_col_sizes(dst_col),& 03114 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 03115 "Unequal column sizes",__LINE__,error) 03116 !CALL dbcsr_unimplemented_error (routineN, "Partial block copy not tested",& 03117 ! error=error, error_level=dbcsr_warning_level) 03118 CALL dbcsr_assert (src_col_offsets(src_col).eq.dst_col_offsets(dst_col),& 03119 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 03120 "Block offsets must be equal",__LINE__,error) 03121 first_col = 1 + MAX (src_col_offsets(src_col), source_start)& 03122 - src_col_offsets(src_col) 03123 last_col = MIN (src_col_offsets(src_col)+ncol-1, source_start+ncol-1) - & 03124 src_col_offsets(src_col+1)-1 + 1 03125 num_col = last_col - first_col + 1 03126 IF (num_col .GT. 0) THEN 03127 CALL dbcsr_block_partial_copy(& 03128 dst_data, row_size, col_size, dst_tr,& 03129 src_data, row_size, col_size, src_tr,& 03130 1, first_col, 1, first_col,& ! offsets 03131 row_size, num_col) ! sizes 03132 ENDIF 03133 CASE (-1) ! Block exists in dst but must be added to dst 03134 IF (source_start .GE. src_col_offsets(src_col) .AND. & 03135 src_col_offsets(src_col)+ncol-1 .LE. source_start+ncol-1) THEN 03136 CALL dbcsr_put_block (matrix_b, src_row, src_col, src_data,& 03137 transposed=src_tr) 03138 ELSE 03139 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_internal_error, routineN,& 03140 "Partial block addition not yet tested",__LINE__,error) 03141 CALL dbcsr_data_clear (tmp_buffer) 03142 CALL dbcsr_block_partial_copy (& 03143 tmp_buffer, row_size, col_size, src_tr,& 03144 src_data, row_size, col_size, src_tr,& 03145 1, first_col, 1, first_col,& 03146 row_size, num_col) 03147 CALL dbcsr_put_block (matrix_b, src_row, src_col, tmp_buffer,& 03148 src_tr) 03149 ENDIF 03150 CASE (1) 03151 ! Ran out of source. Do nothing. 03152 CASE default 03153 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_internal_error,& 03154 routineN, "Trouble syncing iterators",__LINE__,error) 03155 END SELECT 03156 END DO 03157 CALL dbcsr_iterator_stop (src_iter) 03158 CALL dbcsr_iterator_stop (dst_iter) 03159 CALL dbcsr_data_clear_pointer (src_data) 03160 CALL dbcsr_data_clear_pointer (dst_data) 03161 CALL dbcsr_data_release (src_data) 03162 CALL dbcsr_data_release (dst_data) 03163 CALL dbcsr_error_stop(error_handler, error) 03164 END SUBROUTINE dbcsr_copy_columns 03165 03166 03167 ! ***************************************************************************** 03174 ELEMENTAL FUNCTION pos_relation (row1, col1, row2, col2) RESULT (relation) 03175 INTEGER, INTENT(IN) :: row1, col1, row2, col2 03176 INTEGER :: relation 03177 03178 IF (row1 .LT. row2) THEN 03179 relation = -1 03180 ELSEIF (row1 .GT. row2) THEN 03181 relation = 1 03182 ELSE ! rows are equal, check column 03183 IF (col1 .LT. col2) THEN 03184 relation = -1 03185 ELSEIF (col1 .GT. col2) THEN 03186 relation = 1 03187 ELSE 03188 relation = 0 03189 ENDIF 03190 ENDIF 03191 END FUNCTION pos_relation 03192 03193 03194 ! ***************************************************************************** 03208 SUBROUTINE dbcsr_copy_submatrix(matrix_b, matrix_a, name,& 03209 block_row_bounds, block_column_bounds, & 03210 shallow_data, error) 03211 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 03212 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 03213 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name 03214 INTEGER, DIMENSION(2), INTENT(IN), 03215 OPTIONAL :: block_row_bounds, 03216 block_column_bounds 03217 LOGICAL, INTENT(IN), OPTIONAL :: shallow_data 03218 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03219 03220 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_submatrix', 03221 routineP = moduleN//':'//routineN 03222 03223 INTEGER :: blk_p, col, error_handler, 03224 nblocks, new_blk, old_blk, row 03225 INTEGER, ALLOCATABLE, DIMENSION(:) :: blkp_list, col_list, row_list 03226 LOGICAL :: shallow, tr 03227 TYPE(dbcsr_data_obj) :: data_block 03228 TYPE(dbcsr_iterator) :: iter 03229 03230 ! --------------------------------------------------------------------------- 03231 03232 CALL dbcsr_error_set(routineN, error_handler, error) 03233 IF (PRESENT (shallow_data)) THEN 03234 shallow = shallow_data 03235 ELSE 03236 shallow = .FALSE. 03237 ENDIF 03238 CALL dbcsr_access_flush (matrix_a, error=error) 03239 ! Verify assumptions. 03240 IF (PRESENT (block_row_bounds)) THEN 03241 CALL dbcsr_assert (SIZE(block_row_bounds), "EQ", 2, dbcsr_fatal_level,& 03242 dbcsr_wrong_args_error, routineN,& 03243 "Size of bounds specifier must be 2", __LINE__, error=error) 03244 ENDIF 03245 IF (PRESENT (block_column_bounds)) THEN 03246 CALL dbcsr_assert (SIZE(block_column_bounds), "EQ", 2,& 03247 dbcsr_fatal_level,& 03248 dbcsr_wrong_args_error, routineN,& 03249 "Size of bounds specifier must be 2", __LINE__, error=error) 03250 ENDIF 03251 ! Setup target matrix 03252 CALL dbcsr_create (matrix_b, name=name, template=matrix_a, error=error) 03253 CALL dbcsr_finalize(matrix_b, error=error) 03254 IF (.NOT. shallow) THEN 03255 ! Non-shallow copy uses the standard iterator on the source and 03256 ! block put on the target. 03257 ! 03258 !$OMP PARALLEL DEFAULT (none) & 03259 !$OMP PRIVATE (data_block, iter, row, col, tr) & 03260 !$OMP SHARED (matrix_a, matrix_b, error,& 03261 !$OMP block_row_bounds, block_column_bounds) 03262 CALL dbcsr_work_create (matrix_b, work_mutable=.FALSE., error=error) 03263 CALL dbcsr_data_init (data_block) 03264 CALL dbcsr_data_new (data_block, dbcsr_get_data_type (matrix_a)) 03265 CALL dbcsr_iterator_start (iter, matrix_a, dynamic=.TRUE.,& 03266 dynamic_byrows = .TRUE.) 03267 DO WHILE (dbcsr_iterator_blocks_left (iter)) 03268 CALL dbcsr_iterator_next_block (iter, row, col, data_block, tr) 03269 ! Only keep the block if they are within the specified bounds. 03270 IF (PRESENT (block_row_bounds)) THEN 03271 IF (row .LT. block_row_bounds(1)) CYCLE 03272 IF (row .GT. block_row_bounds(2)) CYCLE 03273 ENDIF 03274 IF (PRESENT (block_column_bounds)) THEN 03275 IF (col .LT. block_column_bounds(1)) CYCLE 03276 IF (col .GT. block_column_bounds(2)) CYCLE 03277 ENDIF 03278 CALL dbcsr_put_block (matrix_b, row, col, data_block, transposed=tr) 03279 END DO 03280 CALL dbcsr_iterator_stop (iter) 03281 CALL dbcsr_data_clear_pointer (data_block) 03282 CALL dbcsr_data_release (data_block) 03283 CALL dbcsr_finalize (matrix_b, error=error) 03284 !$OMP END PARALLEL 03285 ELSE 03286 ! For the shallow copy the source matrix data is referenced. 03287 CALL dbcsr_switch_data_area (matrix_b, matrix_a%m%data_area, error=error) 03288 nblocks = dbcsr_get_num_blocks (matrix_a) ! High estimate. 03289 ! Shallow copy goes through source's data blocks and inserts 03290 ! the only the ones corresponding to the submatrix specifier 03291 ! into the target. Block pointers must remain the same as in 03292 ! the source. 03293 ALLOCATE (row_list(nblocks), col_list(nblocks), blkp_list(nblocks)) 03294 ! 03295 CALL dbcsr_iterator_start (iter, matrix_a) 03296 new_blk = 1 03297 DO WHILE (dbcsr_iterator_blocks_left (iter)) 03298 CALL dbcsr_iterator_next_block (iter, row, col,& 03299 blk=old_blk, blk_p=blk_p) 03300 ! Only keep the block if they are within the specified bounds. 03301 IF (PRESENT (block_row_bounds)) THEN 03302 IF (row .LT. block_row_bounds(1)) CYCLE 03303 IF (row .GT. block_row_bounds(2)) CYCLE 03304 ENDIF 03305 IF (PRESENT (block_column_bounds)) THEN 03306 IF (col .LT. block_column_bounds(1)) CYCLE 03307 IF (col .GT. block_column_bounds(2)) CYCLE 03308 ENDIF 03309 row_list(new_blk) = row 03310 col_list(new_blk) = col 03311 blkp_list(new_blk) = blk_p 03312 new_blk = new_blk + 1 03313 ENDDO 03314 new_blk = new_blk - 1 03315 CALL dbcsr_iterator_stop (iter) 03316 CALL dbcsr_reserve_blocks (matrix_b, row_list(1:new_blk),& 03317 col_list(1:new_blk), blkp_list(1:new_blk), error=error) 03318 ENDIF 03319 ! 03320 CALL dbcsr_error_stop(error_handler, error) 03321 END SUBROUTINE dbcsr_copy_submatrix 03322 03323 03324 ! ***************************************************************************** 03337 SUBROUTINE dbcsr_crop_matrix(matrix_b, matrix_a, name,& 03338 full_row_bounds, full_column_bounds, & 03339 shallow_data, error) 03340 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 03341 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 03342 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name 03343 INTEGER, DIMENSION(2), INTENT(IN), 03344 OPTIONAL :: full_row_bounds, 03345 full_column_bounds 03346 LOGICAL, INTENT(IN), OPTIONAL :: shallow_data 03347 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03348 03349 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_crop_matrix', 03350 routineP = moduleN//':'//routineN 03351 03352 INTEGER :: col, error_handler, f_col_f, 03353 f_row_f, l_col_l, l_row_l, row 03354 INTEGER, DIMENSION(2) :: block_col_bounds, 03355 block_row_bounds 03356 LOGICAL :: part_col, part_f_col, 03357 part_f_row, part_l_col, 03358 part_l_row, part_row, 03359 shallow, tr 03360 TYPE(dbcsr_data_obj) :: data_block 03361 TYPE(dbcsr_iterator) :: iter 03362 03363 ! --------------------------------------------------------------------------- 03364 03365 CALL dbcsr_error_set(routineN, error_handler, error) 03366 IF (PRESENT (shallow_data)) THEN 03367 shallow = shallow_data 03368 ELSE 03369 shallow = .FALSE. 03370 ENDIF 03371 CALL dbcsr_access_flush (matrix_a, error=error) 03372 block_row_bounds = 0 03373 block_col_bounds = 0 03374 part_col = .FALSE. 03375 part_row = .FALSE. 03376 ! 03377 ! If row bounds are present, they must be converted to block 03378 ! addressing. 03379 IF (PRESENT (full_row_bounds)) THEN 03380 CALL dbcsr_assert (SIZE(full_row_bounds), "EQ", 2, dbcsr_fatal_level,& 03381 dbcsr_wrong_args_error, routineN,& 03382 "Size of bounds specifier must be 2", __LINE__, error=error) 03383 CALL dbcsr_assert (full_row_bounds(1), "GE", 0,& 03384 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 03385 "Invalid first row bound.", __LINE__, error=error) 03386 CALL dbcsr_assert (full_row_bounds(2), "LE", dbcsr_nfullrows_total(matrix_a),& 03387 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 03388 "Invalid last row bound.", __LINE__, error=error) 03389 part_f_row = .FALSE. 03390 IF (full_row_bounds(1) .EQ. 0) THEN 03391 block_row_bounds(1) = 1 03392 ELSE 03393 CALL find_block_of_element (full_row_bounds(1), block_row_bounds(1),& 03394 dbcsr_nblkrows_total (matrix_a),& 03395 array_data (dbcsr_row_block_offsets (matrix_a)),& 03396 hint = 0, error=error) 03397 part_f_row = array_get (dbcsr_row_block_offsets (matrix_a), block_row_bounds(1))& 03398 .NE. full_row_bounds(1) 03399 ENDIF 03400 f_row_f = -7 03401 IF (part_f_row) THEN 03402 ! Block offset of last cleared row 03403 f_row_f = full_row_bounds(1) -& 03404 array_get (dbcsr_row_block_offsets (matrix_a), block_row_bounds(1)) 03405 ENDIF 03406 part_l_row = .FALSE. 03407 IF (full_row_bounds(2) .EQ. 0) THEN 03408 block_row_bounds(2) = dbcsr_nblkrows_total (matrix_a) 03409 ELSE 03410 CALL find_block_of_element (full_row_bounds(2), block_row_bounds(2),& 03411 dbcsr_nblkrows_total (matrix_a),& 03412 array_data (dbcsr_row_block_offsets (matrix_a)),& 03413 hint = 0, error=error) 03414 part_l_row = array_get (dbcsr_row_block_offsets (matrix_a), block_row_bounds(2)+1)-1& 03415 .NE. full_row_bounds(2) 03416 ENDIF 03417 ! Block offset of first cleared row 03418 l_row_l = -7 03419 IF (part_l_row) THEN 03420 l_row_l = 2 + full_row_bounds(2) - & 03421 array_get (dbcsr_row_block_offsets (matrix_a), block_row_bounds(2)) 03422 ENDIF 03423 part_row = part_f_row .OR. part_l_row 03424 ENDIF 03425 ! 03426 ! If column bounds are present, they must be converted to block 03427 ! addressing. 03428 IF (PRESENT (full_column_bounds)) THEN 03429 CALL dbcsr_assert (SIZE(full_column_bounds), "EQ", 2,& 03430 dbcsr_fatal_level,& 03431 dbcsr_wrong_args_error, routineN,& 03432 "Size of bounds specifier must be 2", __LINE__, error=error) 03433 CALL dbcsr_assert (full_column_bounds(1), "GE", 0,& 03434 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 03435 "Invalid first column bound.", __LINE__, error=error) 03436 CALL dbcsr_assert (full_column_bounds(2), "LE", dbcsr_nfullcols_total(matrix_a),& 03437 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 03438 "Invalid last column bound.", __LINE__, error=error) 03439 part_f_col = .FALSE. 03440 IF (full_column_bounds(1) .EQ. 0) THEN 03441 block_col_bounds(1) = 1 03442 ELSE 03443 CALL find_block_of_element (full_column_bounds(1), block_col_bounds(1),& 03444 dbcsr_nblkcols_total (matrix_a),& 03445 array_data (dbcsr_col_block_offsets (matrix_a)),& 03446 hint=0, error=error) 03447 part_f_col = array_get (dbcsr_col_block_offsets (matrix_a), block_col_bounds(1))& 03448 .NE. full_column_bounds(1) 03449 ENDIF 03450 f_col_f = -7 03451 IF (part_f_col) THEN 03452 ! Block offset of last cleared column 03453 f_col_f = full_column_bounds(1) -& 03454 array_get (dbcsr_col_block_offsets (matrix_a), block_col_bounds(1)) 03455 ENDIF 03456 part_l_col = .FALSE. 03457 IF (full_column_bounds(2) .EQ. 0) THEN 03458 block_col_bounds(2) = dbcsr_nblkcols_total (matrix_a) 03459 ELSE 03460 CALL find_block_of_element (full_column_bounds(2), block_col_bounds(2),& 03461 dbcsr_nblkcols_total (matrix_a),& 03462 array_data (dbcsr_col_block_offsets (matrix_a)),& 03463 hint=0, error=error) 03464 part_l_col = array_get (dbcsr_col_block_offsets (matrix_a), block_col_bounds(2)+1)-1& 03465 .NE. full_column_bounds(2) 03466 ENDIF 03467 l_col_l = -7 03468 IF (part_l_col) THEN 03469 ! Block offset of first cleared column 03470 l_col_l = 2 + full_column_bounds(2) - & 03471 array_get (dbcsr_col_block_offsets (matrix_a), block_col_bounds(2)) 03472 ENDIF 03473 part_col = part_f_col .OR. part_l_col 03474 ENDIF 03475 ! 03476 ! First copy the blocks then perform the intra-block zeroing. 03477 CALL dbcsr_copy_submatrix (matrix_b, matrix_a,& 03478 block_row_bounds=block_row_bounds,& 03479 block_column_bounds=block_col_bounds,& 03480 shallow_data=shallow, error=error) 03481 IF (part_row .OR. part_col) THEN 03482 !$OMP PARALLEL DEFAULT (none) & 03483 !$OMP PRIVATE (data_block, iter, row, col, tr) & 03484 !$OMP SHARED (matrix_b, error,& 03485 !$OMP part_row, part_f_row, part_l_row, f_row_f, l_row_l, & 03486 !$OMP part_col, part_f_col, part_l_col, f_col_f, l_col_l,& 03487 !$OMP block_row_bounds, block_col_bounds) 03488 CALL dbcsr_data_init (data_block) 03489 CALL dbcsr_data_new (data_block, dbcsr_type_1d_to_2d (dbcsr_get_data_type (matrix_b))) 03490 CALL dbcsr_iterator_start (iter, matrix_b, & 03491 dynamic=.TRUE., dynamic_byrows = .TRUE.) 03492 DO WHILE (dbcsr_iterator_blocks_left (iter)) 03493 CALL dbcsr_iterator_next_block (iter, row, col, data_block, tr) 03494 IF (part_row) THEN 03495 IF (row .LT. block_row_bounds(1)) CYCLE 03496 IF (row .GT. block_row_bounds(2)) CYCLE 03497 ENDIF 03498 IF (part_col) THEN 03499 IF (col .LT. block_col_bounds(1)) CYCLE 03500 IF (col .GT. block_col_bounds(2)) CYCLE 03501 ENDIF 03502 IF (part_row) THEN 03503 IF (part_f_row .AND. row .EQ. block_row_bounds(1)) THEN 03504 CALL dbcsr_data_clear (data_block, ub=f_row_f, tr=tr) 03505 ENDIF 03506 IF (part_l_row .AND. row .EQ. block_row_bounds(2)) THEN 03507 CALL dbcsr_data_clear (data_block, lb=l_row_l, tr=tr) 03508 ENDIF 03509 ENDIF 03510 IF (part_col) THEN 03511 IF (part_f_col .AND. col .EQ. block_col_bounds(1)) THEN 03512 CALL dbcsr_data_clear (data_block, ub2=f_col_f, tr=tr) 03513 ENDIF 03514 IF (part_l_col .AND. col .EQ. block_col_bounds(2)) THEN 03515 CALL dbcsr_data_clear (data_block, lb2=l_col_l, tr=tr) 03516 ENDIF 03517 ENDIF 03518 END DO 03519 CALL dbcsr_iterator_stop (iter) 03520 CALL dbcsr_data_clear_pointer (data_block) 03521 CALL dbcsr_data_release (data_block) 03522 CALL dbcsr_finalize (matrix_b, error=error) 03523 !$OMP END PARALLEL 03524 ENDIF 03525 ! 03526 CALL dbcsr_error_stop(error_handler, error) 03527 END SUBROUTINE dbcsr_crop_matrix 03528 03529 03530 03531 03532 ! ***************************************************************************** 03538 SUBROUTINE dbcsr_btriu(matrix_b, matrix_a, error) 03539 03540 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_b 03541 TYPE(dbcsr_obj), INTENT(IN) :: matrix_a 03542 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03543 03544 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_btriu', 03545 routineP = moduleN//':'//routineN 03546 03547 INTEGER :: blk, blk_nze, col, col_size, 03548 error_handler, row, row_size 03549 TYPE(dbcsr_iterator) :: iter 03550 03551 ! --------------------------------------------------------------------------- 03552 03553 CALL dbcsr_error_set(routineN, error_handler, error) 03554 CALL dbcsr_copy(matrix_b, matrix_a, name="triu of "//matrix_a%m%name, error=error) 03555 03556 CALL dbcsr_iterator_start(iter, matrix_b) 03557 03558 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03559 CALL dbcsr_iterator_next_block(iter, row, col, blk=blk,& 03560 row_size=row_size, col_size=col_size) 03561 blk_nze = row_size * col_size 03562 IF(row.GT.col) CALL dbcsr_remove_block(matrix_b, row, col, blk_nze, blk) 03563 ENDDO 03564 03565 CALL dbcsr_iterator_stop(iter) 03566 03567 CALL dbcsr_finalize (matrix_b, error=error) 03568 CALL dbcsr_error_stop(error_handler, error) 03569 END SUBROUTINE dbcsr_btriu 03570 03571 ! ***************************************************************************** 03576 SUBROUTINE dbcsr_triu(matrix_a, error) 03577 03578 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 03579 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03580 03581 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_triu', 03582 routineP = moduleN//':'//routineN 03583 03584 INTEGER :: blk, blk_nze, col, col_size, 03585 error_handler, i, j, row, 03586 row_size 03587 LOGICAL :: tr 03588 REAL(dp), DIMENSION(:, :), POINTER :: DATA 03589 TYPE(dbcsr_iterator) :: iter 03590 03591 ! --------------------------------------------------------------------------- 03592 03593 CALL dbcsr_error_set(routineN, error_handler, error) 03594 CALL dbcsr_iterator_start(iter, matrix_a) 03595 03596 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03597 CALL dbcsr_iterator_next_block(iter, row, col, DATA, tr, & 03598 block_number=blk, row_size=row_size, col_size=col_size) 03599 blk_nze = row_size * col_size 03600 IF(row.GT.col) CALL dbcsr_remove_block(matrix_a, row, col, blk_nze, blk) 03601 IF(row.EQ.col) THEN 03602 DO j = 1,col_size 03603 DO i = j+1,row_size 03604 DATA(i,j) = 0.0_dp 03605 ENDDO 03606 ENDDO 03607 ENDIF 03608 ENDDO 03609 03610 CALL dbcsr_iterator_stop(iter) 03611 03612 CALL dbcsr_finalize (matrix_a, error=error) 03613 CALL dbcsr_error_stop(error_handler, error) 03614 END SUBROUTINE dbcsr_triu 03615 03616 ! ***************************************************************************** 03621 SUBROUTINE dbcsr_symmetrize_block_diag(matrix_a, error) 03622 03623 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 03624 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03625 03626 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_symmetrize_block_diag', 03627 routineP = moduleN//':'//routineN 03628 03629 INTEGER :: blk, col, col_size, 03630 error_handler, i, j, row, 03631 row_size 03632 LOGICAL :: tr 03633 REAL(dp) :: dum 03634 REAL(dp), DIMENSION(:, :), POINTER :: DATA 03635 TYPE(dbcsr_iterator) :: iter 03636 03637 ! --------------------------------------------------------------------------- 03638 03639 CALL dbcsr_error_set(routineN, error_handler, error) 03640 CALL dbcsr_iterator_start(iter, matrix_a) 03641 03642 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03643 CALL dbcsr_iterator_next_block(iter, row, col, DATA, tr, & 03644 block_number=blk, row_size=row_size, col_size=col_size) 03645 IF(col.NE.row) CYCLE 03646 DO j = 1,col_size 03647 DO i = j+1,row_size 03648 dum = (DATA(i,j) + DATA(j,i)) / 2.0_dp 03649 DATA(i,j) = dum 03650 DATA(j,i) = dum 03651 ENDDO 03652 ENDDO 03653 ENDDO 03654 03655 CALL dbcsr_iterator_stop(iter) 03656 03657 CALL dbcsr_finalize(matrix_a, error=error) 03658 CALL dbcsr_error_stop(error_handler, error) 03659 03660 END SUBROUTINE dbcsr_symmetrize_block_diag 03661 03662 ! ***************************************************************************** 03667 SUBROUTINE dbcsr_tril(matrix_a, error) 03668 03669 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 03670 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03671 03672 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_tril', 03673 routineP = moduleN//':'//routineN 03674 03675 INTEGER :: blk, blk_nze, col, col_size, 03676 error_handler, i, j, row, 03677 row_size 03678 LOGICAL :: tr 03679 REAL(dp), DIMENSION(:, :), POINTER :: DATA 03680 TYPE(dbcsr_iterator) :: iter 03681 03682 ! --------------------------------------------------------------------------- 03683 03684 CALL dbcsr_error_set(routineN, error_handler, error) 03685 CALL dbcsr_iterator_start(iter, matrix_a) 03686 03687 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03688 CALL dbcsr_iterator_next_block(iter, row, col, block=DATA, transposed=tr, & 03689 block_number=blk, row_size=row_size, col_size=col_size) 03690 blk_nze = row_size * col_size 03691 IF(row.GT.col) CALL dbcsr_remove_block(matrix_a, row, col, blk_nze, blk) 03692 IF(row.EQ.col) THEN 03693 DO j = 1,col_size 03694 DO i = 1,j-1 03695 DATA(i,j) = 0.0_dp 03696 ENDDO 03697 ENDDO 03698 ENDIF 03699 ENDDO 03700 03701 CALL dbcsr_iterator_stop(iter) 03702 03703 CALL dbcsr_finalize (matrix_a, error=error) 03704 CALL dbcsr_error_stop(error_handler, error) 03705 END SUBROUTINE dbcsr_tril 03706 03707 ! ***************************************************************************** 03717 SUBROUTINE dbcsr_filter_anytype(matrix, eps, method,& 03718 use_absolute, filter_diag, quick, error) 03719 03720 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 03721 TYPE(dbcsr_scalar_type), INTENT(IN) :: eps 03722 INTEGER, INTENT(IN), OPTIONAL :: method 03723 LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag, 03724 quick 03725 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03726 03727 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_filter_anytype', 03728 routineP = moduleN//':'//routineN 03729 03730 COMPLEX(KIND=real_4), DIMENSION(:), 03731 POINTER :: data_c 03732 COMPLEX(KIND=real_8), DIMENSION(:), 03733 POINTER :: data_z 03734 INTEGER :: blk, blk_nze, col, col_size, 03735 error_handler, my_method, 03736 row, row_size 03737 LOGICAL :: gt0, my_filter_diag, tr 03738 REAL(KIND=real_4) :: nrm_s, sdot 03739 REAL(KIND=real_4), DIMENSION(:), POINTER :: data_s 03740 REAL(KIND=real_8) :: ddot, my_absolute, nrm_d 03741 REAL(KIND=real_8), DIMENSION(:), POINTER :: data_d 03742 REAL(KIND=real_8), EXTERNAL :: DZNRM2, SCNRM2 03743 TYPE(dbcsr_iterator) :: iter 03744 03745 ! --------------------------------------------------------------------------- 03746 03747 CALL dbcsr_error_set(routineN, error_handler, error) 03748 my_method = dbcsr_filter_frobenius 03749 IF(PRESENT(method)) my_method = method 03750 my_absolute = 1.0_dp 03751 IF(PRESENT(use_absolute)) my_absolute = dbcsr_maxabs(matrix) 03752 my_filter_diag = .TRUE. 03753 IF(PRESENT(filter_diag)) my_filter_diag = filter_diag 03754 03755 SELECT CASE (eps%data_type) 03756 CASE (dbcsr_type_real_4) 03757 gt0 = eps%r_sp .GT. 0.0_real_4 03758 CASE (dbcsr_type_real_8) 03759 gt0 = eps%r_dp .GT. 0.0_real_8 03760 CASE (dbcsr_type_complex_4) 03761 gt0 = ABS(eps%c_sp) .GT. 0.0_real_4 03762 CASE (dbcsr_type_complex_8) 03763 gt0 = ABS(eps%c_dp) .GT. 0.0_real_8 03764 CASE default 03765 gt0 = .FALSE. 03766 END SELECT 03767 03768 IF(gt0) THEN 03769 03770 CALL dbcsr_iterator_start(iter, matrix, contiguous_pointers=.TRUE.) 03771 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03772 SELECT CASE (dbcsr_get_data_type(matrix)) 03773 03774 CASE (dbcsr_type_real_4) 03775 CALL dbcsr_iterator_next_block(iter, row, col, data_s, tr, blk, & 03776 row_size, col_size) 03777 blk_nze = row_size * col_size 03778 IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks 03779 SELECT CASE(my_method) 03780 CASE(dbcsr_filter_frobenius) 03781 ! 03782 ! Frobenius based 03783 nrm_s = SQRT(SDOT(SIZE(data_s), data_s(1), 1, data_s(1), 1)) 03784 IF(nrm_s.LT.my_absolute*eps%r_sp) & 03785 CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) 03786 CASE DEFAULT 03787 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 03788 routineN,"Only Frobenius based filtering",__LINE__,error) 03789 END SELECT 03790 03791 CASE (dbcsr_type_real_8) 03792 CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk, & 03793 row_size, col_size) 03794 blk_nze = row_size * col_size 03795 IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks 03796 SELECT CASE(my_method) 03797 CASE(dbcsr_filter_frobenius) 03798 ! 03799 ! Frobenius based 03800 nrm_d = SQRT(DDOT(SIZE(data_d), data_d(1), 1, data_d(1), 1)) 03801 IF(nrm_d.LT.my_absolute*eps%r_dp) & 03802 CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) 03803 CASE DEFAULT 03804 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 03805 routineN,"Only Frobenius based filtering",__LINE__,error) 03806 END SELECT 03807 03808 CASE (dbcsr_type_complex_4) 03809 CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk, & 03810 row_size, col_size) 03811 blk_nze = row_size * col_size 03812 IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks 03813 SELECT CASE(my_method) 03814 CASE(dbcsr_filter_frobenius) 03815 ! 03816 ! Frobenius based 03817 nrm_d = SCNRM2(SIZE(data_c), data_c(1), 1) 03818 IF(nrm_d.LT.my_absolute*eps%r_dp) & 03819 CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) 03820 CASE DEFAULT 03821 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 03822 routineN,"Only Frobenius based filtering",__LINE__,error) 03823 END SELECT 03824 03825 CASE (dbcsr_type_complex_8) 03826 CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk, & 03827 row_size, col_size) 03828 blk_nze = row_size * col_size 03829 IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks 03830 SELECT CASE(my_method) 03831 CASE(dbcsr_filter_frobenius) 03832 ! 03833 ! Frobenius based 03834 nrm_d = DZNRM2(SIZE(data_z), data_z(1), 1) 03835 IF(nrm_d.LT.my_absolute*eps%r_dp) & 03836 CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) 03837 CASE DEFAULT 03838 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 03839 routineN,"Only Frobenius based filtering",__LINE__,error) 03840 END SELECT 03841 03842 CASE DEFAULT 03843 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_internal_error, & 03844 routineN,"Wrong data type",__LINE__,error) 03845 END SELECT 03846 ENDDO 03847 CALL dbcsr_iterator_stop(iter) 03848 IF (PRESENT (quick)) THEN 03849 CALL dbcsr_finalize (matrix, reshuffle=.NOT. quick, error=error) 03850 ELSE 03851 CALL dbcsr_finalize (matrix, reshuffle=.TRUE., error=error) 03852 ENDIF 03853 ENDIF 03854 CALL dbcsr_error_stop(error_handler, error) 03855 END SUBROUTINE dbcsr_filter_anytype 03856 03857 ! ***************************************************************************** 03864 SUBROUTINE dbcsr_norm_anytype(matrix, which_norm, norm_scalar, norm_vector, error) 03865 03866 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 03867 INTEGER, INTENT(IN) :: which_norm 03868 REAL(KIND=real_8), INTENT(OUT), OPTIONAL :: norm_scalar 03869 TYPE(dbcsr_data_obj), INTENT(INOUT), 03870 OPTIONAL :: norm_vector 03871 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 03872 03873 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_norm_anytype', 03874 routineP = moduleN//':'//routineN 03875 03876 INTEGER :: blk, col, col_offset, 03877 error_handler, i, j, row, 03878 row_offset 03879 LOGICAL :: tr 03880 TYPE(dbcsr_data_obj) :: data_a 03881 TYPE(dbcsr_iterator) :: iter 03882 03883 ! --------------------------------------------------------------------------- 03884 03885 CALL dbcsr_error_set(routineN, error_handler, error) 03886 03887 SELECT CASE(which_norm) 03888 03889 CASE(dbcsr_norm_frobenius) 03890 03891 IF(PRESENT(norm_scalar)) norm_scalar = dbcsr_frobenius_norm(matrix) 03892 03893 CASE(dbcsr_norm_maxabsnorm) 03894 03895 IF(PRESENT(norm_scalar)) norm_scalar = dbcsr_maxabs(matrix) 03896 03897 CASE(dbcsr_norm_gershgorin) 03898 03899 IF(PRESENT(norm_scalar)) norm_scalar = dbcsr_gershgorin_norm(matrix) 03900 03901 CASE(dbcsr_norm_column) 03902 03903 IF(PRESENT(norm_vector)) THEN 03904 CALL dbcsr_assert (dbcsr_data_get_type(norm_vector), "EQ",& 03905 dbcsr_get_data_type(matrix), dbcsr_fatal_level,& 03906 dbcsr_wrong_args_error, routineN,& 03907 "Mismatched vector/matrix data types", __LINE__, error=error) 03908 IF(dbcsr_has_symmetry(matrix)) THEN 03909 CALL dbcsr_assert (dbcsr_data_get_size(norm_vector), "GE",& 03910 dbcsr_nfullrows_total(matrix), dbcsr_fatal_level,& 03911 dbcsr_wrong_args_error, routineN, "Passed vector too small",& 03912 __LINE__, error=error) 03913 END IF 03914 CALL dbcsr_assert (dbcsr_data_get_size(norm_vector), "GE",& 03915 dbcsr_nfullcols_total(matrix), dbcsr_fatal_level,& 03916 dbcsr_wrong_args_error, routineN, "Passed vector too small",& 03917 __LINE__, error=error) 03918 CALL dbcsr_data_init (data_a) 03919 CALL dbcsr_data_new (data_a, dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix))) 03920 CALL dbcsr_data_clear (norm_vector) 03921 CALL dbcsr_iterator_start(iter, matrix) 03922 DO WHILE (dbcsr_iterator_blocks_left(iter)) 03923 CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr,& 03924 blk, row_offset=row_offset, col_offset=col_offset) 03925 SELECT CASE (dbcsr_get_data_type(matrix)) 03926 CASE (dbcsr_type_real_4) 03927 IF(dbcsr_has_symmetry(matrix) .AND. row.NE.col) THEN 03928 DO j=1,SIZE(data_a%d%r2_sp,2) 03929 DO i=1,SIZE(data_a%d%r2_sp,1) 03930 norm_vector%d%r_sp(col_offset+j-1) & 03931 = norm_vector%d%r_sp(col_offset+j-1) & 03932 + data_a%d%r2_sp(i,j)**2 03933 norm_vector%d%r_sp(row_offset+i-1) & 03934 = norm_vector%d%r_sp(row_offset+i-1) & 03935 + data_a%d%r2_sp(i,j)**2 03936 ENDDO 03937 ENDDO 03938 ELSE 03939 DO j=1,SIZE(data_a%d%r2_sp,2) 03940 DO i=1,SIZE(data_a%d%r2_sp,1) 03941 norm_vector%d%r_sp(col_offset+j-1) & 03942 = norm_vector%d%r_sp(col_offset+j-1) & 03943 + data_a%d%r2_sp(i,j) * data_a%d%r2_sp(i,j) 03944 ENDDO 03945 ENDDO 03946 ENDIF 03947 CASE (dbcsr_type_real_8) 03948 IF(dbcsr_has_symmetry(matrix) .AND. row.NE.col) THEN 03949 DO j=1,SIZE(data_a%d%r2_dp,2) 03950 DO i=1,SIZE(data_a%d%r2_dp,1) 03951 norm_vector%d%r_dp(col_offset+j-1) & 03952 = norm_vector%d%r_dp(col_offset+j-1) & 03953 + data_a%d%r2_dp(i,j)**2 03954 norm_vector%d%r_dp(row_offset+i-1) & 03955 = norm_vector%d%r_dp(row_offset+i-1) & 03956 + data_a%d%r2_dp(i,j)**2 03957 ENDDO 03958 ENDDO 03959 ELSE 03960 DO j=1,SIZE(data_a%d%r2_dp,2) 03961 DO i=1,SIZE(data_a%d%r2_dp,1) 03962 norm_vector%d%r_dp(col_offset+j-1) & 03963 = norm_vector%d%r_dp(col_offset+j-1) & 03964 + data_a%d%r2_dp(i,j) * data_a%d%r2_dp(i,j) 03965 ENDDO 03966 ENDDO 03967 ENDIF 03968 !CASE (dbcsr_type_complex_4) 03969 !CASE (dbcsr_type_complex_8) 03970 CASE DEFAULT 03971 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 03972 routineN,"Only real values",__LINE__,error) 03973 END SELECT 03974 ENDDO 03975 CALL dbcsr_iterator_stop(iter) 03976 CALL dbcsr_data_clear_pointer (data_a) 03977 CALL dbcsr_data_release (data_a) 03978 SELECT CASE (dbcsr_get_data_type(matrix)) 03979 CASE (dbcsr_type_real_4) 03980 CALL mp_sum(norm_vector%d%r_sp,& 03981 dbcsr_mp_group(dbcsr_distribution_mp(matrix%m%dist))) 03982 norm_vector%d%r_sp = SQRT(norm_vector%d%r_sp) 03983 CASE (dbcsr_type_real_8) 03984 CALL mp_sum(norm_vector%d%r_dp,& 03985 dbcsr_mp_group(dbcsr_distribution_mp(matrix%m%dist))) 03986 norm_vector%d%r_dp = SQRT(norm_vector%d%r_dp) 03987 END SELECT 03988 ENDIF 03989 03990 CASE DEFAULT 03991 03992 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN, & 03993 "this norm is NYI",__LINE__,error) 03994 03995 END SELECT 03996 CALL dbcsr_error_stop(error_handler, error) 03997 END SUBROUTINE dbcsr_norm_anytype 03998 03999 SUBROUTINE dbcsr_norm_r8_vec(matrix, which_norm, norm_vector, error) 04000 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04001 INTEGER, INTENT(IN) :: which_norm 04002 REAL(KIND=real_8), DIMENSION(:), 04003 INTENT(OUT), TARGET :: norm_vector 04004 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04005 04006 REAL(KIND=real_8), DIMENSION(:), POINTER :: v_p 04007 TYPE(dbcsr_data_obj) :: norm_vector_a 04008 04009 CALL dbcsr_data_init (norm_vector_a) 04010 CALL dbcsr_data_new (norm_vector_a, dbcsr_type_real_8) 04011 v_p => norm_vector 04012 CALL dbcsr_data_set_pointer (norm_vector_a, v_p) 04013 CALL dbcsr_norm_anytype (matrix, which_norm, norm_vector=norm_vector_a,& 04014 error=error) 04015 CALL dbcsr_data_clear_pointer (norm_vector_a) 04016 CALL dbcsr_data_release (norm_vector_a) 04017 END SUBROUTINE dbcsr_norm_r8_vec 04018 04019 SUBROUTINE dbcsr_norm_r4_scal(matrix, which_norm, norm_scalar, error) 04020 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04021 INTEGER, INTENT(IN) :: which_norm 04022 REAL(KIND=real_4), INTENT(OUT) :: norm_scalar 04023 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04024 04025 REAL(KIND=real_8) :: norm_r8 04026 04027 CALL dbcsr_norm_anytype (matrix, which_norm, norm_scalar=norm_r8,& 04028 error=error) 04029 norm_scalar = REAL(norm_r8, KIND=real_4) 04030 END SUBROUTINE dbcsr_norm_r4_scal 04031 04032 SUBROUTINE dbcsr_norm_r4_vec(matrix, which_norm, norm_vector, error) 04033 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04034 INTEGER, INTENT(IN) :: which_norm 04035 REAL(KIND=real_4), DIMENSION(:), 04036 INTENT(OUT), TARGET :: norm_vector 04037 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04038 04039 REAL(KIND=real_4), DIMENSION(:), POINTER :: v_p 04040 TYPE(dbcsr_data_obj) :: norm_vector_a 04041 04042 CALL dbcsr_data_init (norm_vector_a) 04043 CALL dbcsr_data_new (norm_vector_a, dbcsr_type_real_4) 04044 v_p => norm_vector 04045 CALL dbcsr_data_set_pointer (norm_vector_a, v_p) 04046 CALL dbcsr_norm_anytype (matrix, which_norm, norm_vector=norm_vector_a,& 04047 error=error) 04048 CALL dbcsr_data_clear_pointer (norm_vector_a) 04049 CALL dbcsr_data_release (norm_vector_a) 04050 END SUBROUTINE dbcsr_norm_r4_vec 04051 04052 ! ***************************************************************************** 04058 FUNCTION dbcsr_gershgorin_norm_r8(matrix) RESULT (norm) 04059 04060 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04061 REAL(KIND=real_8) :: norm 04062 04063 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_gershgorin_norm_r8', 04064 routineP = moduleN//':'//routineN 04065 04066 COMPLEX(KIND=real_4), DIMENSION(:, :), 04067 POINTER :: data_c 04068 COMPLEX(KIND=real_8), DIMENSION(:, :), 04069 POINTER :: data_z 04070 INTEGER :: blk, col, col_offset, i, j, 04071 nc, nr, row, row_offset 04072 LOGICAL :: any_sym, tr 04073 REAL(KIND=real_4), DIMENSION(:, :), 04074 POINTER :: data_r 04075 REAL(KIND=real_8), DIMENSION(:, :), 04076 POINTER :: data_d 04077 REAL(real_8), ALLOCATABLE, DIMENSION(:) :: buff_d 04078 TYPE(dbcsr_error_type) :: error 04079 TYPE(dbcsr_iterator) :: iter 04080 04081 ! --------------------------------------------------------------------------- 04082 04083 nr = dbcsr_nfullrows_total(matrix) 04084 nc = dbcsr_nfullcols_total(matrix) 04085 04086 any_sym = dbcsr_get_matrix_type(matrix).EQ.dbcsr_type_symmetric.OR.& 04087 dbcsr_get_matrix_type(matrix).EQ.dbcsr_type_antisymmetric 04088 04089 CALL dbcsr_assert (nr.EQ.nc, dbcsr_fatal_level, dbcsr_wrong_args_error, & 04090 routineN, "not a square matrix",__LINE__,error) 04091 04092 norm = 0.0_dp 04093 ALLOCATE(buff_d(nr)) 04094 buff_d = 0.0_dp 04095 CALL dbcsr_iterator_start(iter, matrix) 04096 DO WHILE (dbcsr_iterator_blocks_left(iter)) 04097 SELECT CASE (dbcsr_get_data_type(matrix)) 04098 CASE (dbcsr_type_real_4) 04099 CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk, & 04100 row_offset=row_offset, col_offset=col_offset) 04101 DO j=1,SIZE(data_r,2) 04102 DO i=1,SIZE(data_r,1) 04103 buff_d(row_offset+i-1) = buff_d(row_offset+i-1) + ABS(data_r(i,j)) 04104 IF(any_sym.AND.row.NE.col) & 04105 buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_r(i,j)) 04106 ENDDO 04107 ENDDO 04108 CASE (dbcsr_type_real_8) 04109 CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk, & 04110 row_offset=row_offset, col_offset=col_offset) 04111 DO j=1,SIZE(data_d,2) 04112 DO i=1,SIZE(data_d,1) 04113 buff_d(row_offset+i-1) = buff_d(row_offset+i-1) + ABS(data_d(i,j)) 04114 IF(any_sym.AND.row.NE.col) & 04115 buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_d(i,j)) 04116 ENDDO 04117 ENDDO 04118 CASE (dbcsr_type_complex_4) 04119 CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk, & 04120 row_offset=row_offset, col_offset=col_offset) 04121 DO j=1,SIZE(data_c,2) 04122 DO i=1,SIZE(data_c,1) 04123 buff_d(row_offset+i-1) = buff_d(row_offset+i-1) + ABS(data_c(i,j)) 04124 IF(any_sym.AND.row.NE.col) & 04125 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04126 routineN,"Only nonsymmetric matrix so far",__LINE__,error) 04127 ! buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_c(i,j)) 04128 ENDDO 04129 ENDDO 04130 CASE (dbcsr_type_complex_8) 04131 CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk, & 04132 row_offset=row_offset, col_offset=col_offset) 04133 DO j=1,SIZE(data_z,2) 04134 DO i=1,SIZE(data_z,1) 04135 buff_d(row_offset+i-1) = buff_d(row_offset+i-1) + ABS(data_z(i,j)) 04136 IF(any_sym.AND.row.NE.col) & 04137 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04138 routineN,"Only nonsymmetric matrix so far",__LINE__,error) 04139 ! buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_z(i,j)) 04140 ENDDO 04141 ENDDO 04142 CASE DEFAULT 04143 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04144 routineN,"Wrong data type",__LINE__,error) 04145 END SELECT 04146 ENDDO 04147 CALL dbcsr_iterator_stop(iter) 04148 CALL mp_sum(buff_d,dbcsr_mp_group(dbcsr_distribution_mp(matrix%m%dist))) 04149 norm = MAXVAL(buff_d) 04150 DEALLOCATE(buff_d) 04151 04152 END FUNCTION dbcsr_gershgorin_norm_r8 04153 04154 ! ***************************************************************************** 04160 FUNCTION dbcsr_maxabs_r8(matrix) RESULT (norm) 04161 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04162 REAL(real_8) :: norm 04163 04164 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_maxabs_r8', 04165 routineP = moduleN//':'//routineN 04166 04167 COMPLEX(KIND=real_4), DIMENSION(:, :), 04168 POINTER :: data_c 04169 COMPLEX(KIND=real_8), DIMENSION(:, :), 04170 POINTER :: data_z 04171 INTEGER :: blk, col, row 04172 LOGICAL :: tr 04173 REAL(KIND=real_4), DIMENSION(:, :), 04174 POINTER :: data_r 04175 REAL(KIND=real_8), DIMENSION(:, :), 04176 POINTER :: data_d 04177 TYPE(dbcsr_error_type) :: error 04178 TYPE(dbcsr_iterator) :: iter 04179 04180 ! --------------------------------------------------------------------------- 04181 04182 norm = 0.0_dp 04183 CALL dbcsr_iterator_start(iter, matrix) 04184 DO WHILE (dbcsr_iterator_blocks_left(iter)) 04185 SELECT CASE (dbcsr_get_data_type(matrix)) 04186 CASE (dbcsr_type_real_4) 04187 CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk) 04188 norm = MAX(norm,REAL(MAXVAL(ABS(data_r)),dp)) 04189 CASE (dbcsr_type_real_8) 04190 CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk) 04191 norm = MAX(norm,MAXVAL(ABS(data_d))) 04192 CASE (dbcsr_type_complex_4) 04193 CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk) 04194 norm = MAX(norm,REAL(MAXVAL(ABS(data_c)),dp)) 04195 CASE (dbcsr_type_complex_8) 04196 CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk) 04197 norm = MAX(norm,MAXVAL(ABS(data_z))) 04198 CASE DEFAULT 04199 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04200 routineN,"Wrong data type",__LINE__,error) 04201 END SELECT 04202 ENDDO 04203 CALL dbcsr_iterator_stop(iter) 04204 04205 !dmp_max, this fixes a bug in g95 04206 ! -> ambigous interface for mp_max 04207 ! in module dbcsr_message_passing and wrongly 04208 ! exported from the module timings 04209 CALL dmp_max(norm, dbcsr_mp_group(dbcsr_distribution_mp(matrix%m%dist))) 04210 04211 END FUNCTION dbcsr_maxabs_r8 04212 04213 ! ***************************************************************************** 04219 FUNCTION dbcsr_frobenius_norm_r8(matrix, local) RESULT (norm) 04220 04221 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix 04222 LOGICAL, INTENT(in), OPTIONAL :: local 04223 REAL(KIND=real_8) :: norm 04224 04225 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_frobenius_norm_r8', 04226 routineP = moduleN//':'//routineN 04227 04228 COMPLEX(KIND=real_4), DIMENSION(:, :), 04229 POINTER :: data_c 04230 COMPLEX(KIND=real_8), DIMENSION(:, :), 04231 POINTER :: data_z 04232 INTEGER :: blk, col, error_handler, row 04233 LOGICAL :: any_sym, my_local, tr 04234 REAL(KIND=real_4), DIMENSION(:, :), 04235 POINTER :: data_r 04236 REAL(KIND=real_8), DIMENSION(:, :), 04237 POINTER :: data_d 04238 REAL(real_8) :: fac 04239 TYPE(dbcsr_error_type) :: error 04240 TYPE(dbcsr_iterator) :: iter 04241 04242 ! --------------------------------------------------------------------------- 04243 04244 CALL dbcsr_error_set(routineN, error_handler, error) 04245 04246 my_local = .FALSE. 04247 IF(PRESENT(local)) my_local = local 04248 04249 any_sym = dbcsr_get_matrix_type(matrix).EQ.dbcsr_type_symmetric.OR.& 04250 dbcsr_get_matrix_type(matrix).EQ.dbcsr_type_antisymmetric 04251 04252 norm = 0.0_dp 04253 CALL dbcsr_iterator_start(iter, matrix) 04254 DO WHILE (dbcsr_iterator_blocks_left(iter)) 04255 SELECT CASE (dbcsr_get_data_type(matrix)) 04256 CASE (dbcsr_type_real_4) 04257 CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk) 04258 fac = 1.0_dp 04259 IF(any_sym.AND.row.NE.col) fac = 2.0_dp 04260 norm = norm + fac * SUM(data_r**2) 04261 CASE (dbcsr_type_real_8) 04262 CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk) 04263 fac = 1.0_dp 04264 IF(any_sym.AND.row.NE.col) fac = 2.0_dp 04265 norm = norm + fac * SUM(data_d**2) 04266 CASE (dbcsr_type_complex_4) 04267 CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk) 04268 fac = 1.0_dp 04269 IF(any_sym.AND.row.NE.col) & 04270 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04271 routineN,"Only nonsymmetric matrix so far",__LINE__,error) 04272 norm = norm + fac * REAL( SUM(CONJG(data_c)*data_c), KIND=real_8) 04273 CASE (dbcsr_type_complex_8) 04274 CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk) 04275 fac = 1.0_dp 04276 IF(any_sym.AND.row.NE.col) & 04277 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04278 routineN,"Only nonsymmetric matrix so far",__LINE__,error) 04279 norm = norm + fac * SUM(CONJG(data_z)*data_z) 04280 CASE DEFAULT 04281 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_unimplemented_error_nr, & 04282 routineN,"Wrong data type",__LINE__,error) 04283 END SELECT 04284 ENDDO 04285 CALL dbcsr_iterator_stop(iter) 04286 IF(.NOT.my_local) CALL mp_sum(norm,dbcsr_mp_group(dbcsr_distribution_mp(matrix%m%dist))) 04287 norm = SQRT(norm) 04288 04289 CALL dbcsr_error_stop(error_handler, error) 04290 04291 END FUNCTION dbcsr_frobenius_norm_r8 04292 04293 04294 ! ***************************************************************************** 04298 SUBROUTINE dbcsr_sum_replicated (matrix, error) 04299 TYPE(dbcsr_obj), INTENT(inout) :: matrix 04300 TYPE(dbcsr_error_type), INTENT(inout) :: error 04301 04302 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_sum_replicated', 04303 routineP = moduleN//':'//routineN 04304 04305 INTEGER :: comm, error_handler, 04306 index_checksum, mynode, 04307 numnodes 04308 INTEGER, ALLOCATABLE, DIMENSION(:) :: all_checksums 04309 TYPE(dbcsr_mp_obj) :: mp 04310 04311 ! --------------------------------------------------------------------------- 04312 04313 CALL dbcsr_error_set(routineN, error_handler, error) 04314 CALL dbcsr_access_flush (matrix, error=error) 04315 mp = dbcsr_distribution_mp (dbcsr_distribution (matrix)) 04316 comm = dbcsr_mp_group (mp) 04317 numnodes = dbcsr_mp_numnodes (mp) 04318 mynode = dbcsr_mp_mynode (mp) 04319 ! 04320 ALLOCATE(all_checksums(numnodes)) 04321 CALL dbcsr_index_checksum (matrix, index_checksum, error) 04322 CALL mp_allgather (index_checksum, all_checksums, comm) 04323 ! 04324 CALL dbcsr_assert (ALL (all_checksums .EQ. index_checksum),& 04325 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 04326 "Replicated matrices do not all have the same index structure.",& 04327 __LINE__, error=error) 04328 ! 04329 SELECT CASE (dbcsr_data_get_type(matrix%m%data_area)) 04330 CASE (dbcsr_type_real_4) 04331 CALL mp_sum (matrix%m%data_area%d%r_sp, comm) 04332 CASE (dbcsr_type_real_8) 04333 CALL mp_sum (matrix%m%data_area%d%r_dp, comm) 04334 CASE (dbcsr_type_complex_4) 04335 CALL mp_sum (matrix%m%data_area%d%c_sp, comm) 04336 CASE (dbcsr_type_complex_8) 04337 CALL mp_sum (matrix%m%data_area%d%c_dp, comm) 04338 CASE default 04339 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_caller_error,& 04340 routineN, "Incorrect data type", __LINE__, error=error) 04341 END SELECT 04342 ! 04343 CALL dbcsr_error_stop(error_handler, error) 04344 END SUBROUTINE dbcsr_sum_replicated 04345 04346 04347 SUBROUTINE dbcsr_trace_a_b_d(matrix_a, matrix_b, trace, trans_a, trans_b, local_sum, error) 04348 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a, matrix_b 04349 REAL(kind=real_8), INTENT(INOUT) :: trace 04350 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: trans_a, trans_b 04351 LOGICAL, INTENT(IN), OPTIONAL :: local_sum 04352 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04353 04354 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_a_b_d', 04355 routineP = moduleN//':'//routineN 04356 04357 INTEGER :: error_handler 04358 REAL(kind=real_4) :: trace_4 04359 04360 CALL dbcsr_error_set(routineN, error_handler, error) 04361 IF( dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & 04362 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8 .OR. & 04363 dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & 04364 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8 .OR. & 04365 dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & 04366 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN 04367 CALL dbcsr_trace_ab_d(matrix_a, matrix_b, trace, trans_a, trans_b, local_sum, error) 04368 ELSEIF(dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & 04369 dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN 04370 trace_4 = 0.0_real_4 04371 CALL dbcsr_trace_ab_s(matrix_a, matrix_b, trace_4, trans_a, trans_b, local_sum, error) 04372 trace = REAL(trace_4,real_8) 04373 ELSE 04374 CALL dbcsr_assert (.FALSE., dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 04375 "Invalid combination of data type, NYI",__LINE__,error) 04376 ENDIF 04377 CALL dbcsr_error_stop(error_handler, error) 04378 END SUBROUTINE dbcsr_trace_a_b_d 04379 04380 04381 04382 SUBROUTINE dbcsr_trace_a_any(matrix_a, trace, error) 04383 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 04384 TYPE(dbcsr_scalar_type), INTENT(INOUT) :: trace 04385 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04386 04387 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_a_any', 04388 routineP = moduleN//':'//routineN 04389 04390 INTEGER :: error_handle 04391 04392 ! --------------------------------------------------------------------------- 04393 04394 CALL dbcsr_error_set(routineN, error_handle, error) 04395 SELECT CASE (dbcsr_scalar_get_type (trace)) 04396 CASE (dbcsr_type_real_4) 04397 CALL dbcsr_trace (matrix_a, trace%r_sp, error) 04398 CASE (dbcsr_type_real_8) 04399 CALL dbcsr_trace (matrix_a, trace%r_dp, error) 04400 CASE (dbcsr_type_complex_4) 04401 CALL dbcsr_trace (matrix_a, trace%c_sp, error) 04402 CASE (dbcsr_type_complex_8) 04403 CALL dbcsr_trace (matrix_a, trace%c_dp, error) 04404 CASE default 04405 CALL dbcsr_assert (.FALSE., dbcsr_failure_level, dbcsr_caller_error,& 04406 routineN, "Invalid data type.",__LINE__,error) 04407 END SELECT 04408 04409 CALL dbcsr_error_stop(error_handle, error) 04410 END SUBROUTINE dbcsr_trace_a_any 04411 04412 04413 ! ***************************************************************************** 04419 FUNCTION dbcsr_block_in_limits(row, col, block_row_limits, block_column_limits) 04420 INTEGER, INTENT(in) :: row, col 04421 INTEGER, DIMENSION(2), INTENT(in), 04422 OPTIONAL :: block_row_limits, 04423 block_column_limits 04424 LOGICAL :: dbcsr_block_in_limits 04425 04426 dbcsr_block_in_limits = .TRUE. 04427 IF (PRESENT (block_row_limits)) THEN 04428 IF (row .LT. block_row_limits(1)) dbcsr_block_in_limits = .FALSE. 04429 IF (row .GT. block_row_limits(2)) dbcsr_block_in_limits = .FALSE. 04430 ENDIF 04431 IF (PRESENT (block_column_limits)) THEN 04432 IF (col .LT. block_column_limits(1)) dbcsr_block_in_limits = .FALSE. 04433 IF (col .GT. block_column_limits(2)) dbcsr_block_in_limits = .FALSE. 04434 ENDIF 04435 END FUNCTION dbcsr_block_in_limits 04436 04437 ! ***************************************************************************** 04449 SUBROUTINE dbcsr_lanczos_extremal_eig(matrix, max_iter, eps, min_eig, max_eig, & 04450 approx_norm_2, converged, error) 04451 04452 TYPE(dbcsr_obj), INTENT(IN) :: matrix 04453 INTEGER, INTENT(IN) :: max_iter 04454 REAL(dp), INTENT(in), OPTIONAL :: eps 04455 REAL(dp), INTENT(out), OPTIONAL :: min_eig, max_eig, 04456 approx_norm_2 04457 LOGICAL, INTENT(OUT), OPTIONAL :: converged 04458 TYPE(dbcsr_error_type), INTENT(inout) :: error 04459 04460 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_lanczos_extremal_eig', 04461 routineP = moduleN//':'//routineN 04462 04463 INTEGER :: error_handler, i, info, 04464 iwork, lwork, n 04465 LOGICAL :: my_converged 04466 REAL(dp) :: alpha, beta, max_eig_old, 04467 min_eig_old, my_eps, 04468 my_max_eig, my_min_eig, 04469 nrm_f, nrm_v 04470 REAL(dp), ALLOCATABLE, DIMENSION(:) :: work 04471 REAL(dp), DIMENSION(max_iter) :: evals 04472 REAL(dp), DIMENSION(max_iter, max_iter) :: T, T0 04473 TYPE(array_i1d_obj) :: col_blk_size, col_dist 04474 TYPE(dbcsr_distribution_obj) :: dist 04475 TYPE(dbcsr_obj) :: f, v, v0 04476 04477 CALL dbcsr_error_set(routineN, error_handler, error) 04478 04479 CALL dbcsr_get_info(matrix,nfullrows_total=n) 04480 04481 CALL create_bl_distribution (col_dist, col_blk_size, 1, & 04482 dbcsr_mp_npcols(dbcsr_distribution_mp(dbcsr_distribution(matrix)))) 04483 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp (dbcsr_distribution(matrix)),& 04484 dbcsr_distribution_row_dist(dbcsr_distribution(matrix)), col_dist) 04485 04486 CALL dbcsr_init(v) 04487 CALL dbcsr_create(v, 'v', dist, dbcsr_type_no_symmetry, matrix%m%row_blk_size,& 04488 col_blk_size, 0, 0, error=error) 04489 CALL dbcsr_finalize(v, error=error) 04490 04491 CALL dbcsr_init(v0) 04492 CALL dbcsr_create(v0, 'v0', dist, dbcsr_type_no_symmetry, matrix%m%row_blk_size,& 04493 col_blk_size, 0, 0, error=error) 04494 CALL dbcsr_finalize(v0, error=error) 04495 04496 CALL dbcsr_init(f) 04497 CALL dbcsr_create(f, 'f', dist, dbcsr_type_no_symmetry, matrix%m%row_blk_size,& 04498 col_blk_size, 0, 0, error=error) 04499 CALL dbcsr_finalize(f, error=error) 04500 04501 CALL dbcsr_distribution_release (dist) 04502 CALL array_release (col_blk_size) 04503 CALL array_release (col_dist) 04504 04505 lwork = 1+2*max_iter+100 04506 ALLOCATE(work(lwork)) 04507 04508 my_eps = 1.0e-1_dp 04509 IF(PRESENT(eps)) my_eps = eps 04510 04511 min_eig_old = 0.0_dp 04512 max_eig_old = 0.0_dp 04513 T(:,:) = 0.0_dp 04514 ! v = rand(n,1) 04515 CALL dbcsr_init_random(v,error=error) 04516 ! v = v / norm(v) 04517 nrm_v = dbcsr_frobenius_norm(v) 04518 CALL dbcsr_scale(v,1.0_dp/nrm_v,error=error) 04519 ! f = A * v 04520 CALL dbcsr_multiply('N','N',1.0_dp,matrix,v,0.0_dp,f,error=error) 04521 ! alpha = f' * v 04522 CALL dbcsr_trace(f,v,alpha,error=error) 04523 ! f = f - alpha * v 04524 CALL dbcsr_add(f,v,1.0_dp,-alpha,error=error) 04525 T(1,1) = alpha 04526 my_min_eig = alpha; my_max_eig = alpha 04527 my_converged = .FALSE. 04528 DO i = 2,max_iter 04529 ! beta = norm(f) 04530 beta = dbcsr_frobenius_norm(f) 04531 ! v0 = v 04532 CALL dbcsr_copy(v0,v,error=error) 04533 ! v = f / beta 04534 CALL dbcsr_add(v,f,0.0_dp,1.0_dp/beta,error=error) 04535 ! f = A * v 04536 CALL dbcsr_multiply('N','N',1.0_dp,matrix,v,0.0_dp,f,error=error) 04537 ! f = f - beta * v0 04538 CALL dbcsr_add(f,v0,1.0_dp,-beta,error=error) 04539 ! alpha = f' * v 04540 CALL dbcsr_trace(f,v,alpha,error=error) 04541 ! f = f - alpha * v 04542 CALL dbcsr_add(f,v,1.0_dp,-alpha,error=error) 04543 T(i ,i-1) = beta 04544 T(i-1,i ) = beta 04545 T(i ,i ) = alpha 04546 ! 04547 max_eig_old = my_max_eig; min_eig_old = my_min_eig 04548 T0(:,:) = T(:,:) 04549 CALL DSYEVD('N','U',i,T0(1,1),max_iter,evals(1),work(1),lwork,iwork,1,info) 04550 CALL dbcsr_assert(info.EQ.0, dbcsr_fatal_level, dbcsr_internal_error, & 04551 routineN, "DSYEVD", __LINE__, error) 04552 my_max_eig = MAXVAL(evals(1:i)); my_min_eig = MINVAL(evals(1:i)) 04553 !write(*,*) routinen//' i',i,'max_eig',my_max_eig,' min_eig',my_min_eig 04554 IF(ABS(my_max_eig-max_eig_old).LT.my_eps.AND.ABS(my_min_eig-min_eig_old).LT.my_eps) THEN 04555 my_converged = .TRUE. 04556 EXIT 04557 ENDIF 04558 ENDDO 04559 04560 IF(PRESENT(approx_norm_2)) THEN 04561 ! norm(f,2) 04562 nrm_f = dbcsr_frobenius_norm(f) 04563 ! norm(T,2) 04564 T0(:,:) = T(:,:) 04565 CALL DSYEVD('N','U',max_iter,T0(1,1),max_iter,evals(1),work(1),lwork,iwork,1,info) 04566 CALL dbcsr_assert(info.EQ.0, dbcsr_fatal_level, dbcsr_internal_error, & 04567 routineN, "DSYEVD", __LINE__, error) 04568 ! norm(T,2)+norm(f,2) 04569 approx_norm_2 = MAXVAL(ABS(evals)) + nrm_f 04570 ENDIF 04571 04572 IF(PRESENT(min_eig)) min_eig = my_min_eig 04573 IF(PRESENT(max_eig)) max_eig = my_max_eig 04574 IF(PRESENT(converged)) converged = my_converged 04575 04576 DEALLOCATE(work) 04577 CALL dbcsr_release(v) 04578 CALL dbcsr_release(v0) 04579 CALL dbcsr_release(f) 04580 04581 CALL dbcsr_error_stop(error_handler, error) 04582 END SUBROUTINE dbcsr_lanczos_extremal_eig 04583 04584 ! ***************************************************************************** 04590 SUBROUTINE dbcsr_init_lib (group, error) 04591 INTEGER, INTENT(IN) :: group 04592 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04593 04594 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_init_lib', 04595 routineP = moduleN//':'//routineN 04596 04597 INTEGER :: error_handle, mynode, ngpus, 04598 numnode 04599 04600 !TODO: problem: init/finalize are called by cp2k_runs AND f77_interface 04601 04602 IF (is_configured) RETURN 04603 04604 CALL dbcsr_error_set(routineN, error_handle, error) 04605 04606 CALL dbcsr_assert (int_1_size, "EQ", 1,& 04607 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 04608 "Incorrect assumption of an 8-bit integer size!",& 04609 __LINE__, error=error) 04610 CALL dbcsr_assert (int_2_size, "EQ", 2,& 04611 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 04612 "Incorrect assumption of a 16-bit integer size!",& 04613 __LINE__, error=error) 04614 CALL dbcsr_assert (int_4_size, "EQ", 4,& 04615 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 04616 "Incorrect assumption of a 32-bit integer size!",& 04617 __LINE__, error=error) 04618 CALL dbcsr_assert (int_8_size, "EQ", 8,& 04619 dbcsr_fatal_level, dbcsr_internal_error, routineN,& 04620 "Incorrect assumption of a 64-bit integer size!",& 04621 __LINE__, error=error) 04622 04623 CALL dbcsr_init_conf (error) 04624 04625 #if defined (__DBCSR_CUDA) 04626 IF (has_cuda) THEN 04627 CALL mp_environ(numnode, mynode, group) 04628 ngpus = dbcsr_cuda_get_n_devices(error) 04629 !$OMP PARALLEL default(none), shared(mynode, ngpus, error) 04630 IF (has_ma) THEN 04631 CALL dbcsr_cuda_init(card_num=ma_set_gpu_affinity(mynode), error=error) 04632 ELSE 04633 CALL dbcsr_cuda_init(card_num=MOD(mynode, ngpus), error=error) 04634 ENDIF 04635 !$OMP END PARALLEL 04636 ENDIF 04637 #endif 04638 04639 !$OMP PARALLEL default(none), shared(error) 04640 CALL dbcsr_mm_cannon_lib_init(error) 04641 !$OMP END PARALLEL 04642 04643 is_configured = .TRUE. 04644 CALL dbcsr_error_stop (error_handle, error) 04645 END SUBROUTINE dbcsr_init_lib 04646 04647 04648 ! ***************************************************************************** 04654 SUBROUTINE dbcsr_finalize_lib (error) 04655 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04656 04657 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_finalize_lib', 04658 routineP = moduleN//':'//routineN 04659 04660 INTEGER :: error_handle 04661 04662 !TODO: problem: init/finalize are called by cp2k_runs AND f77_interface 04663 04664 IF (.NOT. is_configured) RETURN 04665 CALL dbcsr_error_set(routineN, error_handle, error) 04666 04667 !$omp parallel default(none) shared(error) 04668 CALL dbcsr_mm_cannon_lib_finalize(error) 04669 !$omp end parallel 04670 04671 is_configured = .FALSE. 04672 CALL dbcsr_error_stop (error_handle, error) 04673 END SUBROUTINE dbcsr_finalize_lib 04674 04675 04676 04677 ! ***************************************************************************** 04687 SUBROUTINE dbcsr_insert_blocks(matrix_a, matrix_b, error) 04688 TYPE(dbcsr_obj), INTENT(INOUT) :: matrix_a 04689 TYPE(dbcsr_obj), INTENT(IN) :: matrix_b 04690 TYPE(dbcsr_error_type), INTENT(INOUT) :: error 04691 04692 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_insert_blocks', 04693 routineP = moduleN//':'//routineN 04694 04695 INTEGER :: blk, col, data_type_b, 04696 error_handler, nblkrows, 04697 nblks, row 04698 INTEGER, ALLOCATABLE, DIMENSION(:) :: b_row_i 04699 LOGICAL :: tr 04700 TYPE(dbcsr_data_obj) :: data_block 04701 TYPE(dbcsr_iterator) :: iter 04702 04703 ! --------------------------------------------------------------------------- 04704 04705 CALL dbcsr_error_set(routineN, error_handler, error) 04706 ! Checks for validity 04707 CALL dbcsr_assert (dbcsr_valid_index (matrix_a),& 04708 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 04709 "Target matrix A must be valid.", __LINE__, error) 04710 CALL dbcsr_assert (dbcsr_valid_index (matrix_b),& 04711 dbcsr_fatal_level, dbcsr_wrong_args_error, routineN,& 04712 "Source matrix B must be valid.", __LINE__, error) 04713 ! Reserve the blocks to be added 04714 nblks = dbcsr_get_num_blocks (matrix_b) 04715 nblkrows = dbcsr_nblkrows_total (matrix_b) 04716 ALLOCATE (b_row_i(nblks)) 04717 CALL dbcsr_expand_row_index (matrix_b%m%row_p, b_row_i, nblkrows, nblks) 04718 CALL dbcsr_reserve_blocks (matrix_a, b_row_i, matrix_b%m%col_i, error=error) 04719 DEALLOCATE (b_row_i) 04720 ! Prepare data structures 04721 data_type_b = dbcsr_get_data_type (matrix_b) 04722 ! Now add the blocks 04723 CALL dbcsr_data_init (data_block) 04724 CALL dbcsr_data_new (data_block, data_type_b) 04725 CALL dbcsr_iterator_start(iter, matrix_b) 04726 DO WHILE (dbcsr_iterator_blocks_left(iter)) 04727 CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr, blk) 04728 CALL dbcsr_put_block(matrix_a, row, col, data_block, tr,& 04729 summation=.FALSE.) 04730 ENDDO 04731 CALL dbcsr_iterator_stop(iter) 04732 CALL dbcsr_data_clear_pointer (data_block) 04733 CALL dbcsr_data_release (data_block) 04734 ! 04735 CALL dbcsr_finalize (matrix_a, error=error) 04736 CALL dbcsr_error_stop(error_handler, error) 04737 END SUBROUTINE dbcsr_insert_blocks 04738 04739 04740 04741 #include "dbcsr_operations_d_.F" 04742 #include "dbcsr_operations_z_.F" 04743 #include "dbcsr_operations_s_.F" 04744 #include "dbcsr_operations_c_.F" 04745 04746 END MODULE dbcsr_operations
1.7.3