# 1 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" 1 !--------------------------------------------------------------------------------------------------! ! Copyright (C) by the DBCSR developers group - All rights reserved ! ! This file is part of the DBCSR library. ! ! ! ! For information on the license, see the LICENSE file. ! ! For further information please visit https://dbcsr.cp2k.org ! ! SPDX-License-Identifier: GPL-2.0+ ! !--------------------------------------------------------------------------------------------------! MODULE dbcsr_operations !! Higher-level operations on DBCSR matrices. USE dbcsr_array_types, ONLY: array_data, & array_equality, & array_get USE dbcsr_blas_operations, ONLY: set_larnv_seed USE dbcsr_block_access, ONLY: dbcsr_get_block_p, & dbcsr_put_block, & dbcsr_remove_block, & dbcsr_reserve_blocks, & dbcsr_reserve_block2d USE dbcsr_block_operations, ONLY: dbcsr_block_conjg, & dbcsr_block_partial_copy, & dbcsr_block_real_neg, & dbcsr_block_scale, & dbcsr_data_clear USE dbcsr_config, ONLY: default_resize_factor USE dbcsr_data_methods, ONLY: & dbcsr_data_clear_pointer, dbcsr_data_ensure_size, dbcsr_data_get_size, & dbcsr_data_get_size_referenced, dbcsr_data_get_type, dbcsr_data_init, dbcsr_data_new, & dbcsr_data_release, dbcsr_data_set_pointer, dbcsr_get_data, dbcsr_scalar, & dbcsr_scalar_are_equal, dbcsr_scalar_fill_all, dbcsr_scalar_get_type, dbcsr_scalar_one, & dbcsr_scalar_zero, dbcsr_type_1d_to_2d, dbcsr_get_data_p_d, dbcsr_get_data_p_s, dbcsr_scalar_set_type USE dbcsr_data_operations, ONLY: dbcsr_data_convert, & dbcsr_data_copyall, & dbcsr_switch_data_area USE dbcsr_dist_methods, ONLY: dbcsr_distribution_col_dist, & dbcsr_distribution_local_cols, & dbcsr_distribution_local_rows, & dbcsr_distribution_mp, & dbcsr_distribution_row_dist, & dbcsr_distribution_get USE dbcsr_dist_operations, ONLY: checker_square_proc, & checker_tr, & dbcsr_find_column, & dbcsr_get_stored_coordinates, & dbcsr_get_stored_block_info USE dbcsr_index_operations, ONLY: dbcsr_index_checksum, & dbcsr_index_compact, & dbcsr_repoint_index USE dbcsr_iterator_operations, ONLY: dbcsr_iterator_blocks_left, & dbcsr_iterator_next_block, & dbcsr_iterator_start, & dbcsr_iterator_stop USE dbcsr_methods, ONLY: & dbcsr_col_block_offsets, dbcsr_distribution, dbcsr_get_data_size, dbcsr_get_data_type, & dbcsr_get_index_memory_type, dbcsr_get_matrix_type, dbcsr_get_num_blocks, & dbcsr_get_replication_type, dbcsr_has_symmetry, dbcsr_max_col_size, dbcsr_max_row_size, & dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_nfullcols_total, dbcsr_nfullrows_total, & dbcsr_row_block_offsets, dbcsr_valid_index, dbcsr_get_nze, dbcsr_nfullcols_local, & dbcsr_nfullrows_local, dbcsr_get_num_blocks, dbcsr_release, dbcsr_wm_use_mutable USE dbcsr_mpiwrap, ONLY: mp_comm_type USE dbcsr_mp_methods, ONLY: dbcsr_mp_group, & dbcsr_mp_mynode, & dbcsr_mp_mypcol, & dbcsr_mp_myprow, & dbcsr_mp_numnodes, & dbcsr_mp_pgrid USE dbcsr_ptr_util, ONLY: ensure_array_size, & pointer_view USE dbcsr_toollib, ONLY: swap USE dbcsr_types, ONLY: & dbcsr_data_obj, dbcsr_distribution_obj, dbcsr_filter_frobenius, dbcsr_func_artanh, & dbcsr_func_asin, dbcsr_func_cos, dbcsr_func_ddsin, dbcsr_func_ddtanh, dbcsr_func_dsin, & dbcsr_func_dtanh, dbcsr_func_inverse, dbcsr_func_inverse_special, dbcsr_func_sin, & dbcsr_func_spread_from_zero, dbcsr_func_tanh, dbcsr_func_truncate, dbcsr_iterator, & dbcsr_mp_obj, dbcsr_norm_column, dbcsr_norm_frobenius, dbcsr_norm_gershgorin, & dbcsr_norm_maxabsnorm, dbcsr_repl_full, dbcsr_repl_none, dbcsr_scalar_type, dbcsr_type, & dbcsr_type_antihermitian, dbcsr_type_antisymmetric, dbcsr_type_complex_4, & dbcsr_type_complex_8, dbcsr_type_hermitian, dbcsr_type_no_symmetry, dbcsr_type_real_4, & dbcsr_type_real_8, dbcsr_type_symmetric USE dbcsr_dist_util, ONLY: find_block_of_element USE dbcsr_work_operations, ONLY: dbcsr_create, & dbcsr_finalize, & dbcsr_work_create, & add_work_coordinate USE dbcsr_kinds, ONLY: dp, & int_4, & int_8, & real_4, & real_8, & sp USE dbcsr_mpiwrap, ONLY: mp_allgather, & mp_max, & mp_sum #include "base/dbcsr_base_uses.f90" !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads IMPLICIT NONE PRIVATE CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbcsr_operations' ! prettify protection CHARACTER, PARAMETER :: xa = dbcsr_type_hermitian, xb = dbcsr_type_antihermitian, & xc = dbcsr_type_no_symmetry PUBLIC :: dbcsr_trace, dbcsr_dot, dbcsr_add_on_diag, & dbcsr_set, dbcsr_scale, dbcsr_add, dbcsr_copy, & dbcsr_copy_into_existing, & dbcsr_get_diag, dbcsr_set_diag, & dbcsr_get_block_diag, dbcsr_hadamard_product, & dbcsr_filter, dbcsr_filter_anytype, dbcsr_scale_by_vector, & dbcsr_function_of_elements, & dbcsr_triu, & dbcsr_init_random PUBLIC :: dbcsr_sum_replicated PUBLIC :: dbcsr_norm_scalar, dbcsr_norm_r8_vec, dbcsr_conjg, & dbcsr_gershgorin_norm, dbcsr_maxabs, dbcsr_frobenius_norm PUBLIC :: dbcsr_crop_matrix PUBLIC :: dbcsr_get_info, dbcsr_may_be_dense, dbcsr_get_occupation PUBLIC :: dbcsr_clear, dbcsr_add_block_node, dbcsr_conform_scalar PUBLIC :: dbcsr_zero ! The interfaces for the generic routines found in the generated ! generic files. INTERFACE dbcsr_conform_scalar MODULE PROCEDURE make_conformant_scalar_d, make_conformant_scalar_s MODULE PROCEDURE make_conformant_scalar_c, make_conformant_scalar_z END INTERFACE INTERFACE dbcsr_trace MODULE PROCEDURE dbcsr_trace_s, dbcsr_trace_sd, & dbcsr_trace_c, dbcsr_trace_z END INTERFACE INTERFACE dbcsr_dot MODULE PROCEDURE dbcsr_dot_s, dbcsr_dot_sd, & dbcsr_dot_c, dbcsr_dot_z END INTERFACE INTERFACE dbcsr_scale MODULE PROCEDURE dbcsr_scale_anytype MODULE PROCEDURE dbcsr_scale_s, dbcsr_scale_d, dbcsr_scale_c, dbcsr_scale_z END INTERFACE INTERFACE dbcsr_scale_by_vector MODULE PROCEDURE dbcsr_scale_by_vector_anytype MODULE PROCEDURE dbcsr_scale_by_vector_s, dbcsr_scale_by_vector_d MODULE PROCEDURE dbcsr_scale_by_vector_c, dbcsr_scale_by_vector_z END INTERFACE INTERFACE dbcsr_set MODULE PROCEDURE dbcsr_set_s, dbcsr_set_d, dbcsr_set_c, dbcsr_set_z END INTERFACE INTERFACE dbcsr_add MODULE PROCEDURE dbcsr_add_anytype MODULE PROCEDURE dbcsr_add_s, dbcsr_add_d, dbcsr_add_c, dbcsr_add_z END INTERFACE INTERFACE dbcsr_add_on_diag MODULE PROCEDURE dbcsr_add_on_diag_s, dbcsr_add_on_diag_ds MODULE PROCEDURE dbcsr_add_on_diag_c, dbcsr_add_on_diag_z END INTERFACE INTERFACE dbcsr_filter MODULE PROCEDURE dbcsr_filter_anytype MODULE PROCEDURE dbcsr_filter_s, dbcsr_filter_d, & dbcsr_filter_c, dbcsr_filter_z END INTERFACE INTERFACE dbcsr_get_diag MODULE PROCEDURE dbcsr_get_diag_s, dbcsr_get_diag_d, dbcsr_get_diag_c, dbcsr_get_diag_z END INTERFACE INTERFACE dbcsr_set_diag MODULE PROCEDURE dbcsr_set_diag_s, dbcsr_set_diag_d, dbcsr_set_diag_c, dbcsr_set_diag_z END INTERFACE LOGICAL, PARAMETER :: debug_mod = .FALSE. LOGICAL, PARAMETER :: careful_mod = .FALSE. INTEGER, PARAMETER, PRIVATE :: rpslot_owner = 1 INTEGER, PARAMETER, PRIVATE :: rpslot_addblks = 2 INTEGER, PARAMETER, PRIVATE :: rpslot_addoffset = 3 INTEGER, PARAMETER, PRIVATE :: rpslot_oldblks = 4 INTEGER, PARAMETER, PRIVATE :: rpslot_oldoffset = 5 INTEGER, PARAMETER, PRIVATE :: rpslot_totaloffset = 6 INTEGER, PARAMETER, PRIVATE :: rpnslots = 6 CONTAINS # 1 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" 1 # 9 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 11 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 169 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 197 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" 2 # 198 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" FUNCTION make_conformant_scalar_d (scalar, matrix) RESULT(encapsulated) !! Encapsulates a given scalar value and makes it conform with the !! type of the matrix. REAL(kind=real_8), INTENT(IN) :: scalar TYPE(dbcsr_type), INTENT(IN) :: matrix TYPE(dbcsr_scalar_type) :: encapsulated INTEGER :: data_type, scalar_data_type encapsulated = dbcsr_scalar(scalar) CALL dbcsr_scalar_fill_all(encapsulated) data_type = dbcsr_get_data_type(matrix) scalar_data_type = dbcsr_scalar_get_type(encapsulated) IF (scalar_data_type .EQ. dbcsr_type_complex_4 .OR. & scalar_data_type .EQ. dbcsr_type_complex_8) THEN IF (data_type .NE. dbcsr_type_complex_4 .AND. data_type .NE. dbcsr_type_complex_8) & DBCSR_ABORT("Can not conform a complex to a real number") END IF CALL dbcsr_scalar_set_type(encapsulated, data_type) END FUNCTION make_conformant_scalar_d # 198 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" FUNCTION make_conformant_scalar_s (scalar, matrix) RESULT(encapsulated) !! Encapsulates a given scalar value and makes it conform with the !! type of the matrix. REAL(kind=real_4), INTENT(IN) :: scalar TYPE(dbcsr_type), INTENT(IN) :: matrix TYPE(dbcsr_scalar_type) :: encapsulated INTEGER :: data_type, scalar_data_type encapsulated = dbcsr_scalar(scalar) CALL dbcsr_scalar_fill_all(encapsulated) data_type = dbcsr_get_data_type(matrix) scalar_data_type = dbcsr_scalar_get_type(encapsulated) IF (scalar_data_type .EQ. dbcsr_type_complex_4 .OR. & scalar_data_type .EQ. dbcsr_type_complex_8) THEN IF (data_type .NE. dbcsr_type_complex_4 .AND. data_type .NE. dbcsr_type_complex_8) & DBCSR_ABORT("Can not conform a complex to a real number") END IF CALL dbcsr_scalar_set_type(encapsulated, data_type) END FUNCTION make_conformant_scalar_s # 198 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" FUNCTION make_conformant_scalar_z (scalar, matrix) RESULT(encapsulated) !! Encapsulates a given scalar value and makes it conform with the !! type of the matrix. COMPLEX(kind=real_8), INTENT(IN) :: scalar TYPE(dbcsr_type), INTENT(IN) :: matrix TYPE(dbcsr_scalar_type) :: encapsulated INTEGER :: data_type, scalar_data_type encapsulated = dbcsr_scalar(scalar) CALL dbcsr_scalar_fill_all(encapsulated) data_type = dbcsr_get_data_type(matrix) scalar_data_type = dbcsr_scalar_get_type(encapsulated) IF (scalar_data_type .EQ. dbcsr_type_complex_4 .OR. & scalar_data_type .EQ. dbcsr_type_complex_8) THEN IF (data_type .NE. dbcsr_type_complex_4 .AND. data_type .NE. dbcsr_type_complex_8) & DBCSR_ABORT("Can not conform a complex to a real number") END IF CALL dbcsr_scalar_set_type(encapsulated, data_type) END FUNCTION make_conformant_scalar_z # 198 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" FUNCTION make_conformant_scalar_c (scalar, matrix) RESULT(encapsulated) !! Encapsulates a given scalar value and makes it conform with the !! type of the matrix. COMPLEX(kind=real_4), INTENT(IN) :: scalar TYPE(dbcsr_type), INTENT(IN) :: matrix TYPE(dbcsr_scalar_type) :: encapsulated INTEGER :: data_type, scalar_data_type encapsulated = dbcsr_scalar(scalar) CALL dbcsr_scalar_fill_all(encapsulated) data_type = dbcsr_get_data_type(matrix) scalar_data_type = dbcsr_scalar_get_type(encapsulated) IF (scalar_data_type .EQ. dbcsr_type_complex_4 .OR. & scalar_data_type .EQ. dbcsr_type_complex_8) THEN IF (data_type .NE. dbcsr_type_complex_4 .AND. data_type .NE. dbcsr_type_complex_8) & DBCSR_ABORT("Can not conform a complex to a real number") END IF CALL dbcsr_scalar_set_type(encapsulated, data_type) END FUNCTION make_conformant_scalar_c # 220 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" SUBROUTINE dbcsr_add_block_node(matrix, block_row, block_col, block) !! Emulation of sparse_matrix_types/add_block_node mapped !! to add_real_matrix_block.... should not be used any longer !! It adds a block to the dbcsr matrix and returns a rank-2 pointer to the !! block. Currently it only and always uses the mutable data. TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix INTEGER, INTENT(IN) :: block_row, block_col !! the row !! the column REAL(KIND=dp), DIMENSION(:, :), POINTER :: block !! the block to put INTEGER :: c, ithread, mynode, p, r LOGICAL :: dbg, existed, is_there, tr TYPE(dbcsr_distribution_obj) :: dist ! --------------------------------------------------------------------------- dbg = .FALSE. ithread = 0 !$ ithread = omp_get_thread_num() IF (.NOT. ASSOCIATED(matrix%wms)) THEN CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) matrix%valid = .FALSE. END IF !$ IF (SIZE(matrix%wms) .LT. omp_get_num_threads()) & !$ DBCSR_ABORT("Too few threads.") IF (.NOT. dbcsr_wm_use_mutable(matrix%wms(ithread + 1))) & DBCSR_ABORT("Data loss due to no conversion of appendable to mutable data") is_there = ASSOCIATED(block) !r = row ; c = col ; tr = .FALSE. !CALL dbcsr_get_stored_coordinates (matrix, r, c, tr) !CALL dbcsr_reserve_block2d (matrix, row, col, block) !write(*,*) 'add_block_node: block_row',block_row,' block_col',block_col CALL dbcsr_reserve_block2d(matrix, block_row, block_col, block, & existed=existed) ! IF (dbg) THEN r = block_row; c = block_col; tr = .FALSE. CALL dbcsr_get_stored_coordinates(matrix, r, c, p) CALL dbcsr_get_info(matrix, distribution=dist) CALL dbcsr_distribution_get(dist, mynode=mynode) IF (p .NE. mynode) & DBCSR_WARN("Adding non-local element") END IF IF (existed) DBCSR_WARN("You should not add existing blocks according to old API.") IF (.NOT. is_there) block(:, :) = 0.0_dp END SUBROUTINE dbcsr_add_block_node SUBROUTINE dbcsr_conjg(matrix) !! Conjugate a DBCSR matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_conjg' INTEGER :: blk, col, data_type, handle, row LOGICAL :: tr TYPE(dbcsr_data_obj) :: data_any TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- ! CALL timeset(routineN, handle) data_type = dbcsr_get_data_type(matrix) CALL dbcsr_data_init(data_any) CALL dbcsr_data_new(data_any, data_type) CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, data_any, tr, blk) SELECT CASE (data_type) CASE (dbcsr_type_complex_4) data_any%d%c_sp = CONJG(data_any%d%c_sp) CASE (dbcsr_type_complex_8) data_any%d%c_dp = CONJG(data_any%d%c_dp) CASE DEFAULT ! needed for g95 END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_any) CALL dbcsr_data_release(data_any) CALL timestop(handle) END SUBROUTINE dbcsr_conjg SUBROUTINE dbcsr_zero(matrix_a) !! fill a dbcsr matrix with zeros TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_zero' INTEGER :: handle CALL timeset(routineN, handle) SELECT CASE (dbcsr_get_data_type(matrix_a)) #if defined(__DBCSR_DISABLE_WORKSHARE) CASE (dbcsr_type_complex_4) matrix_a%data_area%d%c_sp = (0.0, 0.0) CASE (dbcsr_type_complex_8) matrix_a%data_area%d%c_dp = (0.0_dp, 0.0_dp) CASE (dbcsr_type_real_4) matrix_a%data_area%d%r_sp = 0.0 CASE (dbcsr_type_real_8) matrix_a%data_area%d%r_dp = 0.0_dp #else CASE (dbcsr_type_complex_4) !$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(matrix_a) matrix_a%data_area%d%c_sp = (0.0, 0.0) !$OMP END PARALLEL WORKSHARE CASE (dbcsr_type_complex_8) !$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(matrix_a) matrix_a%data_area%d%c_dp = (0.0_dp, 0.0_dp) !$OMP END PARALLEL WORKSHARE CASE (dbcsr_type_real_4) !$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(matrix_a) matrix_a%data_area%d%r_sp = 0.0 !$OMP END PARALLEL WORKSHARE CASE (dbcsr_type_real_8) !$OMP PARALLEL WORKSHARE DEFAULT(NONE), SHARED(matrix_a) matrix_a%data_area%d%r_dp = 0.0_dp !$OMP END PARALLEL WORKSHARE #endif END SELECT CALL timestop(handle) END SUBROUTINE dbcsr_zero SUBROUTINE dbcsr_scale_anytype(matrix_a, alpha_scalar, limits) !! Scales a DBCSR matrix by alpha !! !! Limits !! A 4-tuple describing (first_row, last_row, first_column, last_column). Set !! to 0 to avoid limiting. TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: alpha_scalar !! a scalar INTEGER, DIMENSION(4), INTENT(IN), OPTIONAL :: limits !! Scale only a subbox CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_anytype' INTEGER, PARAMETER :: first_col_i = 3, first_row_i = 1, & last_col_i = 4, last_row_i = 2 INTEGER :: a_col, a_col_size, a_row, a_row_size, col_offset, handle, row_offset, & scale_col_offset, scale_col_size, scale_row_offset, scale_row_size INTEGER, DIMENSION(4) :: my_limits LOGICAL :: do_scale, has_limits, tr TYPE(dbcsr_data_obj) :: data_any TYPE(dbcsr_iterator) :: iter TYPE(dbcsr_scalar_type) :: one ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) ! Limits are only honored if the argument is present and any are ! non-zero. IF (PRESENT(limits)) THEN has_limits = ANY(limits(:) .NE. 0) ELSE has_limits = .FALSE. END IF my_limits(first_row_i) = 1 my_limits(last_row_i) = dbcsr_nfullrows_total(matrix_a) my_limits(first_col_i) = 1 my_limits(last_col_i) = dbcsr_nfullcols_total(matrix_a) IF (has_limits) THEN IF (limits(last_col_i) .NE. 0) THEN IF (debug_mod .AND. (limits(last_col_i) < 0 .OR. limits(last_col_i) > dbcsr_nfullcols_total(matrix_a))) & DBCSR_ABORT("Specified last column is out of bounds.") my_limits(last_col_i) = limits(last_col_i) END IF IF (limits(first_col_i) .NE. 0) THEN IF (debug_mod .AND. (limits(first_col_i) < 0 .OR. limits(first_col_i) > dbcsr_nfullcols_total(matrix_a))) & DBCSR_ABORT("Specified first column is out of bounds.") my_limits(first_col_i) = limits(first_col_i) END IF IF (limits(last_row_i) .NE. 0) THEN IF (debug_mod .AND. (limits(last_row_i) < 0 .OR. limits(last_row_i) > dbcsr_nfullrows_total(matrix_a))) & DBCSR_ABORT("Specified last row is out of bounds.") my_limits(last_row_i) = limits(last_row_i) END IF IF (limits(first_row_i) .NE. 0) THEN IF (debug_mod .AND. (limits(first_row_i) < 0 .OR. limits(first_row_i) > dbcsr_nfullrows_total(matrix_a))) & DBCSR_ABORT("Specified first row is out of bounds.") my_limits(first_row_i) = limits(first_row_i) END IF END IF ! ! quick return if possible one = dbcsr_scalar_one(dbcsr_scalar_get_type(alpha_scalar)) do_scale = .NOT. dbcsr_scalar_are_equal(alpha_scalar, one) ! ! let's go IF (do_scale) THEN !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (iter, data_any) & !$OMP PRIVATE (a_row, a_col, tr, a_row_size, a_col_size, & !$OMP row_offset, col_offset) & !$OMP PRIVATE (scale_row_size, scale_col_size,& !$OMP scale_row_offset, scale_col_offset) & !$OMP SHARED (matrix_a, my_limits,alpha_scalar) CALL dbcsr_data_init(data_any) CALL dbcsr_data_new(data_any, dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix_a))) CALL dbcsr_iterator_start(iter, matrix_a, read_only=.FALSE., & contiguous_pointers=.FALSE., dynamic=.TRUE., & dynamic_byrows=.TRUE., shared=.TRUE.) iterations: DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, a_row, a_col, data_any, tr, & row_size=a_row_size, col_size=a_col_size, & row_offset=row_offset, col_offset=col_offset) IF (a_row_size .GT. 0 .AND. a_col_size .GT. 0) THEN CALL frame_block_limit(a_row_size, row_offset, & my_limits(first_row_i), my_limits(last_row_i), & scale_row_size, scale_row_offset) CALL frame_block_limit(a_col_size, col_offset, & my_limits(first_col_i), my_limits(last_col_i), & scale_col_size, scale_col_offset) IF (tr) THEN CALL swap(scale_row_size, scale_col_size) CALL swap(scale_row_offset, scale_col_offset) END IF CALL dbcsr_block_scale(data_any, scale=alpha_scalar, & row_size=scale_row_size, col_size=scale_col_size, & lb=scale_row_offset, lb2=scale_col_offset) END IF END DO iterations CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_any) CALL dbcsr_data_release(data_any) !$OMP END PARALLEL END IF CALL timestop(handle) END SUBROUTINE dbcsr_scale_anytype ELEMENTAL SUBROUTINE frame_block_limit(block_size, block_offset, & first_limit, last_limit, & frame_size, frame_offset) !! Determines the effect of limits on a block INTEGER, INTENT(IN) :: block_size, block_offset, first_limit, & last_limit !! size of block !! global offset of block !! lower limit !! upper limit INTEGER, INTENT(OUT) :: frame_size, frame_offset !! size of block region within the limits !! starting position of the block region that is within the limits INTEGER :: f, l f = MAX(block_offset, first_limit) l = MIN(block_offset + block_size - 1, last_limit) frame_size = MAX(l - f + 1, 0) frame_offset = MIN(f - block_offset + 1, block_size) END SUBROUTINE frame_block_limit SUBROUTINE dbcsr_scale_by_vector_anytype(matrix_a, alpha, side) !! Scales a DBCSR matrix by alpha TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_data_obj), INTENT(IN), OPTIONAL :: alpha !! the scaling vector CHARACTER(LEN=*), INTENT(IN) :: side !! apply the scaling from the side CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_by_vector_anytype' INTEGER :: a_blk, a_col, a_col_size, a_nze, a_row, & a_row_size, col_offset, data_type, & handle, i, icol, irow, row_offset LOGICAL :: right, tr TYPE(dbcsr_data_obj) :: data_any TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) ! check that alpha and matrix have the same data type IF (dbcsr_get_data_type(matrix_a) .NE. alpha%d%data_type) & DBCSR_ABORT("wrong data type matrix_a") IF (ASSOCIATED(alpha%d%r2_sp) .OR. ASSOCIATED(alpha%d%r2_dp) & .OR. ASSOCIATED(alpha%d%c2_sp) .OR. ASSOCIATED(alpha%d%c2_dp)) & DBCSR_ABORT("alpha is not a vector") ! ! set vars right = .TRUE. SELECT CASE (side) CASE ('right'); right = .TRUE. CASE ('left'); right = .FALSE. CASE DEFAULT DBCSR_ABORT("wrong side="//side) END SELECT ! check that alpha and matrix have matching sizes IF (right .AND. dbcsr_nfullcols_total(matrix_a) /= dbcsr_data_get_size(alpha)) THEN DBCSR_ABORT("vector size does not match matrix row size for RHS scaling") ELSE IF ((.NOT. right) .AND. dbcsr_nfullrows_total(matrix_a) /= dbcsr_data_get_size(alpha)) THEN DBCSR_ABORT("vector size does not match matrix col size for LHS scaling") END IF ! ! let's go data_type = dbcsr_get_data_type(matrix_a) CALL dbcsr_data_init(data_any) CALL dbcsr_data_new(data_any, dbcsr_get_data_type(matrix_a)) CALL dbcsr_iterator_start(iter, matrix_a) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, a_row, a_col, data_any, tr, & block_number=a_blk, & row_size=a_row_size, col_size=a_col_size, & row_offset=row_offset, col_offset=col_offset) a_nze = a_row_size*a_col_size IF (a_nze .EQ. 0) CYCLE ! Skip empty blocks ! ! let's scale IF (right) THEN SELECT CASE (data_type) CASE (dbcsr_type_real_4) DO i = 1, a_col_size DO icol = (i - 1)*a_row_size + 1, (i - 1)*a_row_size + a_row_size data_any%d%r_sp(icol) = data_any%d%r_sp(icol)*alpha%d%r_sp(col_offset + i - 1) END DO END DO CASE (dbcsr_type_real_8) DO i = 1, a_col_size DO icol = (i - 1)*a_row_size + 1, (i - 1)*a_row_size + a_row_size data_any%d%r_dp(icol) = data_any%d%r_dp(icol)*alpha%d%r_dp(col_offset + i - 1) END DO END DO CASE (dbcsr_type_complex_4) DO i = 1, a_col_size DO icol = (i - 1)*a_row_size + 1, (i - 1)*a_row_size + a_row_size data_any%d%c_sp(icol) = data_any%d%c_sp(icol)*alpha%d%c_sp(col_offset + i - 1) END DO END DO CASE (dbcsr_type_complex_8) DO i = 1, a_col_size DO icol = (i - 1)*a_row_size + 1, (i - 1)*a_row_size + a_row_size data_any%d%c_dp(icol) = data_any%d%c_dp(icol)*alpha%d%c_dp(col_offset + i - 1) END DO END DO END SELECT ELSE SELECT CASE (data_type) CASE (dbcsr_type_real_4) DO i = 1, a_row_size DO irow = i, i + a_col_size*a_row_size - 1, a_row_size data_any%d%r_sp(irow) = data_any%d%r_sp(irow)*alpha%d%r_sp(row_offset + i - 1) END DO END DO CASE (dbcsr_type_real_8) DO i = 1, a_row_size DO irow = i, i + a_col_size*a_row_size - 1, a_row_size data_any%d%r_dp(irow) = data_any%d%r_dp(irow)*alpha%d%r_dp(row_offset + i - 1) END DO END DO CASE (dbcsr_type_complex_4) DO i = 1, a_row_size DO irow = i, i + a_col_size*a_row_size - 1, a_row_size data_any%d%c_sp(irow) = data_any%d%c_sp(irow)*alpha%d%c_sp(row_offset + i - 1) END DO END DO CASE (dbcsr_type_complex_8) DO i = 1, a_row_size DO irow = i, i + a_col_size*a_row_size - 1, a_row_size data_any%d%c_dp(irow) = data_any%d%c_dp(irow)*alpha%d%c_dp(row_offset + i - 1) END DO END DO END SELECT END IF END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_any) CALL dbcsr_data_release(data_any) CALL timestop(handle) END SUBROUTINE dbcsr_scale_by_vector_anytype SUBROUTINE dbcsr_add_anytype(matrix_a, matrix_b, alpha_scalar, beta_scalar, flop) !! add and scale matrices !! A = alpha*A + beta*B or TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN), OPTIONAL :: alpha_scalar, beta_scalar INTEGER(KIND=int_8), INTENT(INOUT), OPTIONAL :: flop CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_anytype' INTEGER :: data_type_a, data_type_b, & handle, size_work, iw INTEGER(KIND=int_8) :: my_flop, local_matrix_size LOGICAL :: do_scale TYPE(dbcsr_iterator) :: iter TYPE(dbcsr_scalar_type) :: my_beta_scalar ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (.NOT. dbcsr_valid_index(matrix_a)) & DBCSR_ABORT("Invalid matrix") IF ((dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_symmetric .OR. & dbcsr_get_matrix_type(matrix_a) .EQ. dbcsr_type_antisymmetric) .NEQV. & (dbcsr_get_matrix_type(matrix_b) .EQ. dbcsr_type_symmetric .OR. & dbcsr_get_matrix_type(matrix_b) .EQ. dbcsr_type_antisymmetric)) THEN DBCSR_ABORT("Summing general with symmetric matrix NYI") END IF data_type_a = dbcsr_get_data_type(matrix_a) data_type_b = dbcsr_get_data_type(matrix_b) ! my_beta_scalar = dbcsr_scalar_one(data_type_b) IF (PRESENT(beta_scalar)) my_beta_scalar = beta_scalar ! ! let's go IF ((dbcsr_nblkrows_total(matrix_a) .NE. dbcsr_nblkrows_total(matrix_b)) .OR. & (dbcsr_nblkcols_total(matrix_a) .NE. dbcsr_nblkcols_total(matrix_b)) .OR. & (data_type_a .NE. data_type_b)) & DBCSR_ABORT("matrices not consistent") IF (data_type_a .NE. my_beta_scalar%data_type) & DBCSR_ABORT("beta type parameter not consistent with matrices type") do_scale = .NOT. dbcsr_scalar_are_equal(my_beta_scalar, dbcsr_scalar_one(data_type_b)) IF (PRESENT(alpha_scalar)) THEN CALL dbcsr_scale(matrix_a, alpha_scalar=alpha_scalar) END IF IF ((.NOT. dbcsr_scalar_are_equal(my_beta_scalar, & dbcsr_scalar_zero(data_type_b))) .AND. & dbcsr_get_num_blocks(matrix_b) .GT. 0) THEN ! Pre-size work arrays of matrix_a to avoid continuous reallocation. ! Overestimate for symmetric matrix and multiple threads! local_matrix_size = INT(dbcsr_nfullrows_local(matrix_a), KIND=int_8)* & dbcsr_nfullcols_local(matrix_a) size_work = MAX(0, INT(MIN(local_matrix_size - INT(dbcsr_get_nze(matrix_a), KIND=int_8), & INT(dbcsr_get_nze(matrix_b), KIND=int_8)), KIND=int_4)) my_flop = 0 !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (iter, iw) & !$OMP SHARED (matrix_a, matrix_b, data_type_b, size_work) & !$OMP SHARED (do_scale, my_beta_scalar) & !$OMP REDUCTION (+ : my_flop) CALL dbcsr_work_create(matrix_a, & nblks_guess=matrix_b%nblks, & sizedata_guess=size_work, & work_mutable=.FALSE.) !$OMP BARRIER iw = 1 !$ iw = omp_get_thread_num() + 1 CALL dbcsr_iterator_start(iter, matrix_b, & shared=.TRUE., read_only=.TRUE., contiguous_pointers=.FALSE., & dynamic=.TRUE., dynamic_byrows=.TRUE.) SELECT CASE (data_type_b) CASE (dbcsr_type_real_4) CALL dbcsr_add_anytype_s(matrix_a, matrix_b, iter, iw, do_scale, my_beta_scalar, my_flop) CASE (dbcsr_type_real_8) CALL dbcsr_add_anytype_d(matrix_a, matrix_b, iter, iw, do_scale, my_beta_scalar, my_flop) CASE (dbcsr_type_complex_4) CALL dbcsr_add_anytype_c(matrix_a, matrix_b, iter, iw, do_scale, my_beta_scalar, my_flop) CASE (dbcsr_type_complex_8) CALL dbcsr_add_anytype_z(matrix_a, matrix_b, iter, iw, do_scale, my_beta_scalar, my_flop) CASE default DBCSR_ABORT("Invalid data type") END SELECT CALL dbcsr_iterator_stop(iter) CALL dbcsr_finalize(matrix_a) !$OMP END PARALLEL IF (PRESENT(flop)) flop = flop + my_flop END IF CALL timestop(handle) END SUBROUTINE dbcsr_add_anytype SUBROUTINE dbcsr_add_d(matrix_a, matrix_b, alpha_scalar, beta_scalar) !! Interface for dbcsr_add TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a TYPE(dbcsr_type), INTENT(IN) :: matrix_b REAL(real_8), INTENT(IN) :: alpha_scalar, beta_scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_d' INTEGER :: handle CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(alpha_scalar), & beta_scalar=dbcsr_scalar(beta_scalar)) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(REAL(alpha_scalar, real_4)), & beta_scalar=dbcsr_scalar(REAL(beta_scalar, real_4))) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_4) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(CMPLX(alpha_scalar, 0, real_4)), & beta_scalar=dbcsr_scalar(CMPLX(beta_scalar, 0, real_4))) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_8 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_8) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(CMPLX(alpha_scalar, 0, real_8)), & beta_scalar=dbcsr_scalar(CMPLX(beta_scalar, 0, real_8))) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_add_d SUBROUTINE dbcsr_add_s(matrix_a, matrix_b, alpha_scalar, beta_scalar) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a TYPE(dbcsr_type), INTENT(IN) :: matrix_b REAL(real_4), INTENT(IN) :: alpha_scalar, beta_scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_s' INTEGER :: handle CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(alpha_scalar), & beta_scalar=dbcsr_scalar(beta_scalar)) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_add_s SUBROUTINE dbcsr_add_z(matrix_a, matrix_b, alpha_scalar, beta_scalar) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a TYPE(dbcsr_type), INTENT(IN) :: matrix_b COMPLEX(real_8), INTENT(IN) :: alpha_scalar, beta_scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_z' INTEGER :: handle CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_8 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_8) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(alpha_scalar), & beta_scalar=dbcsr_scalar(beta_scalar)) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_4) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(CMPLX(alpha_scalar, KIND=real_4)), & beta_scalar=dbcsr_scalar(CMPLX(beta_scalar, KIND=real_4))) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_add_z SUBROUTINE dbcsr_add_c(matrix_a, matrix_b, alpha_scalar, beta_scalar) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a TYPE(dbcsr_type), INTENT(IN) :: matrix_b COMPLEX(real_4), INTENT(IN) :: alpha_scalar, beta_scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_c' INTEGER :: handle CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_complex_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_complex_4) THEN CALL dbcsr_add_anytype(matrix_a, matrix_b, & alpha_scalar=dbcsr_scalar(alpha_scalar), & beta_scalar=dbcsr_scalar(beta_scalar)) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_add_c SUBROUTINE dbcsr_add_on_diag_ds(matrix, alpha) TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_8), INTENT(IN) :: alpha IF (dbcsr_get_data_type(matrix) == dbcsr_type_real_4) THEN CALL dbcsr_add_on_diag_s(matrix, REAL(alpha, kind=real_4)) ELSE CALL dbcsr_add_on_diag_d(matrix, alpha) END IF END SUBROUTINE dbcsr_add_on_diag_ds SUBROUTINE dbcsr_function_of_elements(matrix_a, func, a0, a1, a2) !! Computes various functions (defined by func) of matrix elements !! @note sign(A,B) returns the value of A with the sign of B !! dbcsr_func_inverse: 1/(a1*x+a0) !! fails if the inversion produces infinite numbers !! dbcsr_func_inverse_special: 1/(x+sign(a0,x)) !! safe inverse: if a0>0 then the denominator is never zero !! dbcsr_func_tanh: tanh(a1*x+a0) !! dbcsr_func_dtanh: d(tanh(a1*x+a0)) / dx !! dbcsr_func_ddtanh: d2(tanh(a1*x+a0)) / dx2 !! dbcsr_func_artanh: artanh(a1*x+a0)=ln[(1+(a1*x+a0))/(1-(a1*x+a0))]/2 !! fails if |a1*x+a0| >= 1 !! dbcsr_func_sread_from_zero: if |x|<|a0| then x=sign(a0,x) !! dbcsr_func_truncate: if |x|>|a0| then x=sign(a0,x) !! dbcsr_func_sin: sin(a1*x+a0) !! dbcsr_func_cos: cos(a1*x+a0) !! dbcsr_func_dsin: d(sin(a1*x+a0)) / dx = a1*cos(a1*x+a0) !! dbcsr_func_ddsin: d2(sin(a1*x+a0)) / dx2 = -a1*a1*sin(a1*x+a0) !! dbcsr_func_asin: asin(a1*x+a0) !! fails if |a1*x+a0| > 1 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix INTEGER, INTENT(IN) :: func REAL(kind=dp), INTENT(IN), OPTIONAL :: a0, a1, a2 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_function_of_elements' INTEGER :: blk, col, col_size, data_type, handle, & ielem, nze, row, row_size LOGICAL :: tr_a REAL(kind=dp) :: p0, p1, p2 TYPE(dbcsr_data_obj) :: a_data TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (PRESENT(a0)) THEN p0 = a0 ELSE p0 = 0.0_dp END IF IF (PRESENT(a1)) THEN p1 = a1 ELSE p1 = 1.0_dp END IF IF (PRESENT(a2)) THEN p2 = a2 ELSE p2 = 0.0_dp END IF data_type = dbcsr_get_data_type(matrix_a) CALL dbcsr_data_init(a_data) CALL dbcsr_data_new(a_data, data_type) CALL dbcsr_iterator_start(iter, matrix_a) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, a_data, tr_a, blk, & row_size=row_size, col_size=col_size) nze = row_size*col_size SELECT CASE (data_type) !CASE (dbcsr_type_real_4) ! a_data%d%r_sp(1:nze) = 1.0_real_4/a_data%d%r_sp(1:nze) ! IF(MAXVAL(ABS(a_data%d%r_sp)).GE.HUGE(0.0_real_4))& ! DBCSR_ABORT("Division by zero") CASE (dbcsr_type_real_8) SELECT CASE (func) CASE (dbcsr_func_spread_from_zero) ! if |x|<|a0| then x=|a0|*sign(x) DO ielem = 1, nze IF (ABS(a_data%d%r_dp(ielem)) .LT. ABS(p0)) THEN a_data%d%r_dp(ielem) = SIGN(p0, a_data%d%r_dp(ielem)) END IF END DO CASE (dbcsr_func_truncate) ! if |x|>|a0| then x=|a0|*sign(x) DO ielem = 1, nze IF (ABS(a_data%d%r_dp(ielem)) .GT. ABS(p0)) THEN a_data%d%r_dp(ielem) = SIGN(p0, a_data%d%r_dp(ielem)) END IF END DO CASE (dbcsr_func_inverse_special) !IF (MINVAL(ABS(a_data%d%r_dp)).le.ABS(p2)) THEN ! ! there is at least one near-zero element, ! ! invert element-by-element ! DO ielem=1,nze ! IF (a_data%d%r_dp(ielem).le.ABS(p2)) THEN ! a_data%d%r_dp(ielem) = 0.0_real_8 ! ELSE ! a_data%d%r_dp(ielem) = & ! 1.0_real_8/(p1*a_data%d%r_dp(ielem)+p0) ! ENDIF ! ENDDO !ELSE ! a_data%d%r_dp(1:nze) = 1.0_real_8/(p1*a_data%d%r_dp(1:nze)+p0) !ENDIF a_data%d%r_dp(1:nze) = 1.0_real_8/(a_data%d%r_dp(1:nze) + SIGN(p0, a_data%d%r_dp(1:nze))) CASE (dbcsr_func_inverse) a_data%d%r_dp(1:nze) = 1.0_real_8/(p1*a_data%d%r_dp(1:nze) + p0) IF (MAXVAL(ABS(a_data%d%r_dp)) .GE. HUGE(0.0_real_8)) & DBCSR_ABORT("Division by zero") CASE (dbcsr_func_tanh) a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0) CASE (dbcsr_func_dtanh) a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0) a_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)**2 a_data%d%r_dp(1:nze) = p1*(1.0_real_8 - a_data%d%r_dp(1:nze)) CASE (dbcsr_func_ddtanh) a_data%d%r_dp(1:nze) = TANH(p1*a_data%d%r_dp(1:nze) + p0) a_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)**3 - a_data%d%r_dp(1:nze) a_data%d%r_dp(1:nze) = 2.0_real_8*(p1**2)*a_data%d%r_dp(1:nze) CASE (dbcsr_func_artanh) a_data%d%r_dp(1:nze) = p1*a_data%d%r_dp(1:nze) + p0 IF (MAXVAL(ABS(a_data%d%r_dp)) .GE. 1.0_real_8) & DBCSR_ABORT("ARTANH is undefined for |x|>=1") a_data%d%r_dp(1:nze) = (1.0_real_8 + a_data%d%r_dp(1:nze)) & /(1.0_real_8 - a_data%d%r_dp(1:nze)) a_data%d%r_dp(1:nze) = 0.5_real_8*LOG(a_data%d%r_dp(1:nze)) CASE (dbcsr_func_sin) a_data%d%r_dp(1:nze) = SIN(p1*a_data%d%r_dp(1:nze) + p0) CASE (dbcsr_func_cos) a_data%d%r_dp(1:nze) = COS(p1*a_data%d%r_dp(1:nze) + p0) CASE (dbcsr_func_dsin) a_data%d%r_dp(1:nze) = p1*COS(p1*a_data%d%r_dp(1:nze) + p0) CASE (dbcsr_func_ddsin) a_data%d%r_dp(1:nze) = -p1*p1*SIN(p1*a_data%d%r_dp(1:nze) + p0) CASE (dbcsr_func_asin) a_data%d%r_dp(1:nze) = p1*a_data%d%r_dp(1:nze) + p0 IF (MAXVAL(ABS(a_data%d%r_dp)) .GT. 1.0_real_8) & DBCSR_ABORT("ASIN is undefined for |x|>1") a_data%d%r_dp(1:nze) = ASIN(a_data%d%r_dp(1:nze)) CASE DEFAULT DBCSR_ABORT("Unknown function of matrix elements") END SELECT !CASE (dbcsr_type_complex_4) !CASE (dbcsr_type_complex_8) CASE DEFAULT DBCSR_ABORT("Operation is implemented only for dp real values") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(a_data) CALL dbcsr_data_release(a_data) CALL timestop(handle) END SUBROUTINE dbcsr_function_of_elements SUBROUTINE dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c, & b_assume_value) !! Hadamard product !! C = A . B (C needs to be different from A and B) TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrix !! DBCSR matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix_c !! DBCSR matrix REAL(KIND=dp), INTENT(IN), OPTIONAL :: b_assume_value CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_hadamard_product' INTEGER :: blk, col, col_size, data_type, handle, & nze, row, row_size LOGICAL :: assume_blocks_in_b, found, tr_a, tr_b REAL(KIND=dp) :: assumed_b_value TYPE(dbcsr_data_obj) :: a_data, b_data, c_data TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- IF (PRESENT(b_assume_value)) THEN assume_blocks_in_b = .TRUE. assumed_b_value = b_assume_value ELSE assume_blocks_in_b = .FALSE. assumed_b_value = 0.0_dp END IF CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .NE. dbcsr_get_data_type(matrix_b) .OR. & dbcsr_get_data_type(matrix_a) .NE. dbcsr_get_data_type(matrix_c)) & DBCSR_ABORT("data types not consistent, need to fix that") IF (dbcsr_nblkrows_total(matrix_a) .NE. dbcsr_nblkrows_total(matrix_b) .OR. & dbcsr_nblkrows_total(matrix_c) .NE. dbcsr_nblkrows_total(matrix_a)) & DBCSR_ABORT("matrices not consistent") data_type = dbcsr_get_data_type(matrix_a) CALL dbcsr_data_init(c_data) CALL dbcsr_data_new(c_data, data_type, & data_size=dbcsr_max_row_size(matrix_a)*dbcsr_max_col_size(matrix_a)) CALL dbcsr_zero(matrix_c) CALL dbcsr_data_init(a_data) CALL dbcsr_data_new(a_data, data_type) CALL dbcsr_data_init(b_data) CALL dbcsr_data_new(b_data, data_type) CALL dbcsr_iterator_start(iter, matrix_a) DO WHILE (dbcsr_iterator_blocks_left(iter)) SELECT CASE (dbcsr_get_data_type(matrix_a)) !CASE (dbcsr_type_real_4) CASE (dbcsr_type_real_8) CALL dbcsr_iterator_next_block(iter, row, col, a_data, tr_a, blk, & row_size=row_size, col_size=col_size) nze = row_size*col_size CALL dbcsr_get_block_p(matrix_b, row, col, b_data, tr_b, found) IF (tr_a .NEQV. tr_b) & DBCSR_ABORT("tr not consistent, need to fix that") IF (found) THEN SELECT CASE (data_type) CASE (dbcsr_type_real_4) c_data%d%r_sp(1:nze) = a_data%d%r_sp(1:nze)*b_data%d%r_sp(1:nze) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%r_sp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_real_8) c_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)*b_data%d%r_dp(1:nze) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%r_dp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_complex_4) c_data%d%c_sp(1:nze) = a_data%d%c_sp(1:nze)*b_data%d%c_sp(1:nze) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%c_sp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_complex_8) c_data%d%c_dp(1:nze) = a_data%d%c_dp(1:nze)*b_data%d%c_dp(1:nze) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%c_dp(1:nze), transposed=tr_a, & summation=.FALSE.) END SELECT ELSE IF (assume_blocks_in_b) THEN ! this makes not too much sense, to delete ? SELECT CASE (data_type) CASE (dbcsr_type_real_4) c_data%d%r_sp(1:nze) = a_data%d%r_sp(1:nze)*REAL(assumed_b_value, KIND=sp) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%r_sp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_real_8) c_data%d%r_dp(1:nze) = a_data%d%r_dp(1:nze)*assumed_b_value CALL dbcsr_put_block(matrix_c, row, col, c_data%d%r_dp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_complex_4) c_data%d%c_sp(1:nze) = a_data%d%c_sp(1:nze)*REAL(assumed_b_value, KIND=sp) CALL dbcsr_put_block(matrix_c, row, col, c_data%d%c_sp(1:nze), transposed=tr_a, & summation=.FALSE.) CASE (dbcsr_type_complex_8) c_data%d%c_dp(1:nze) = a_data%d%c_dp(1:nze)*assumed_b_value CALL dbcsr_put_block(matrix_c, row, col, c_data%d%c_dp(1:nze), transposed=tr_a, & summation=.FALSE.) END SELECT END IF END IF !CASE (dbcsr_type_complex_4) !CASE (dbcsr_type_complex_8) CASE DEFAULT DBCSR_ABORT("Only real double precision") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_finalize(matrix_c) CALL dbcsr_data_clear_pointer(a_data) CALL dbcsr_data_clear_pointer(b_data) CALL dbcsr_data_release(c_data) CALL dbcsr_data_release(a_data) CALL dbcsr_data_release(b_data) CALL timestop(handle) END SUBROUTINE dbcsr_hadamard_product SUBROUTINE dbcsr_init_random(matrix, keep_sparsity, mini_seed) !! ... TODO : unify with other version which is generic in the data_type TYPE(dbcsr_type), INTENT(INOUT) :: matrix LOGICAL, OPTIONAL :: keep_sparsity INTEGER, INTENT(IN), OPTIONAL :: mini_seed CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_init_random' INTEGER :: col, col_size, handle, hold, iseed(4), & mynode, ncol, nrow, row, row_size, & stored_col, stored_row, my_mini_seed INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size LOGICAL :: found, my_keep_sparsity, tr REAL(real_8), ALLOCATABLE, DIMENSION(:) :: rnd REAL(real_8), DIMENSION(:, :), POINTER :: buff, data_d ! --------------------------------------------------------------------------- my_keep_sparsity = .FALSE. IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity my_mini_seed = 1 IF (PRESENT(mini_seed)) my_mini_seed = mini_seed CALL timeset(routineN, handle) row_blk_size => array_data(matrix%row_blk_size) col_blk_size => array_data(matrix%col_blk_size) mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) ALLOCATE (rnd(MAXVAL(row_blk_size)*MAXVAL(col_blk_size))) nrow = dbcsr_nblkrows_total(matrix) ncol = dbcsr_nblkcols_total(matrix) DO row = 1, nrow DO col = 1, ncol row_size = row_blk_size(row) col_size = col_blk_size(col) tr = .FALSE. stored_row = row stored_col = col CALL dbcsr_get_stored_coordinates(matrix, stored_row, stored_col, hold) IF (hold .EQ. mynode) THEN CALL dbcsr_get_block_p(matrix, stored_row, stored_col, data_d, tr, found) IF (found .OR. (.NOT. my_keep_sparsity)) THEN ! set the seed for dlarnv, is here to guarantee same value of the random numbers ! for all layouts (and block distributions) CALL set_larnv_seed(row, nrow, col, ncol, my_mini_seed, iseed) CALL dlarnv(1, iseed, row_size*col_size, rnd(1)) END IF IF (found) THEN CALL dcopy(row_size*col_size, rnd, 1, data_d, 1) ELSE IF (.NOT. my_keep_sparsity) THEN ALLOCATE (buff(row_size, col_size)) CALL dcopy(row_size*col_size, rnd, 1, buff, 1) CALL dbcsr_put_block(matrix, stored_row, stored_col, buff) DEALLOCATE (buff) END IF END IF END IF END DO END DO DEALLOCATE (rnd) CALL dbcsr_finalize(matrix) CALL timestop(handle) END SUBROUTINE dbcsr_init_random SUBROUTINE dbcsr_get_block_diag(matrix, diag) !! get the diagonal of a dbcsr matrix TYPE(dbcsr_type), INTENT(IN) :: matrix !! the matrix TYPE(dbcsr_type), INTENT(INOUT) :: diag !! the diagonal CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_block_diag' INTEGER :: blk, col, handle, row LOGICAL :: tr TYPE(dbcsr_data_obj) :: data_a TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) CALL dbcsr_create(diag, name='diag of '//TRIM(matrix%name), & template=matrix) CALL dbcsr_data_init(data_a) CALL dbcsr_data_new(data_a, dbcsr_get_data_type(matrix)) CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr, blk) IF (row .EQ. col) CALL dbcsr_put_block(diag, row, col, data_a, transposed=tr) END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_a) CALL dbcsr_data_release(data_a) CALL dbcsr_finalize(diag) CALL timestop(handle) END SUBROUTINE dbcsr_get_block_diag LOGICAL FUNCTION symmetry_consistent(matrix_type, data_type) !! checks if matrix symmetry and data_type are consistent !! \brief note: does not check the symmetry of the data itself CHARACTER, INTENT(IN) :: matrix_type INTEGER, INTENT(IN) :: data_type symmetry_consistent = .FALSE. SELECT CASE (data_type) CASE (dbcsr_type_real_4, dbcsr_type_real_8) SELECT CASE (matrix_type) CASE (dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_type_antisymmetric) symmetry_consistent = .TRUE. END SELECT CASE (dbcsr_type_complex_4, dbcsr_type_complex_8) SELECT CASE (matrix_type) CASE (dbcsr_type_no_symmetry, dbcsr_type_hermitian, dbcsr_type_antihermitian) symmetry_consistent = .TRUE. END SELECT CASE DEFAULT DBCSR_ABORT("Invalid data type.") END SELECT END FUNCTION symmetry_consistent LOGICAL FUNCTION symmetry_compatible(matrix_type_a, matrix_type_b) !! checks if symmetries of two matrices are compatible for copying !! \brief data from matrix_a(source) to matrix_b(target) CHARACTER, INTENT(IN) :: matrix_type_a, matrix_type_b symmetry_compatible = .FALSE. SELECT CASE (matrix_type_a) CASE (dbcsr_type_no_symmetry) SELECT CASE (matrix_type_b) CASE (dbcsr_type_no_symmetry) symmetry_compatible = .TRUE. END SELECT CASE (dbcsr_type_symmetric, dbcsr_type_hermitian) SELECT CASE (matrix_type_b) CASE (dbcsr_type_symmetric, dbcsr_type_hermitian) symmetry_compatible = .TRUE. END SELECT CASE (dbcsr_type_antisymmetric, dbcsr_type_antihermitian) SELECT CASE (matrix_type_b) CASE (dbcsr_type_antisymmetric, dbcsr_type_antihermitian) symmetry_compatible = .TRUE. END SELECT CASE DEFAULT DBCSR_ABORT("Invalid matrix type.") END SELECT END FUNCTION symmetry_compatible SUBROUTINE dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, & shallow_data, keep_imaginary, matrix_type) !! copy a matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b !! target DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! source DBCSR matrix CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name !! name of the new matrix LOGICAL, INTENT(IN), OPTIONAL :: keep_sparsity, shallow_data, & keep_imaginary !! keep the target matrix sparsity; default is False. !! shallow data copy !! when copy from complex to real,& the default is to keep only the real part; if this flag is set, the imaginary part is !! used CHARACTER, INTENT(IN), OPTIONAL :: matrix_type !! 'N' for normal, 'T' for transposed, 'S' for symmetric, and 'A' for antisymmetric CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy' CHARACTER :: new_matrix_type, repl_type INTEGER :: handle, new_type LOGICAL :: keep_sparse, shallow ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (.NOT. symmetry_consistent(dbcsr_get_matrix_type(matrix_a), dbcsr_get_data_type(matrix_a))) & DBCSR_ABORT("Source matrix symmetry not consistent with its data type.") shallow = .FALSE.; IF (PRESENT(shallow_data)) shallow = shallow_data keep_sparse = .FALSE. IF (PRESENT(keep_sparsity)) keep_sparse = keep_sparsity IF (keep_sparse .AND. .NOT. dbcsr_valid_index(matrix_b)) & DBCSR_ABORT("Target matrix must be valid to keep its sparsity") IF (keep_sparse .AND. shallow) & DBCSR_WARN("Shallow copy not compatibly with sparsity retainment") IF (keep_sparse) THEN IF (PRESENT(name)) matrix_b%name = name CALL dbcsr_copy_into_existing(matrix_b, matrix_a) ELSE IF (dbcsr_valid_index(matrix_b)) THEN new_type = dbcsr_get_data_type(matrix_b) repl_type = dbcsr_get_replication_type(matrix_b) ELSE new_type = dbcsr_get_data_type(matrix_a) repl_type = dbcsr_get_replication_type(matrix_a) END IF new_matrix_type = dbcsr_get_matrix_type(matrix_a) IF (PRESENT(matrix_type)) THEN IF (.NOT. symmetry_compatible(dbcsr_get_matrix_type(matrix_a), matrix_type)) & CALL dbcsr_abort(__LOCATION__, "Specified target matrix symmetry "//matrix_type// & " not compatible with source matrix type "//dbcsr_get_matrix_type(matrix_a)) new_matrix_type = matrix_type END IF IF (.NOT. symmetry_consistent(new_matrix_type, new_type)) & CALL dbcsr_abort(__LOCATION__, "Target matrix symmetry "// & new_matrix_type//" not consistent with its data type.") IF (PRESENT(name)) THEN CALL dbcsr_create(matrix_b, name=TRIM(name), & template=matrix_a, & matrix_type=new_matrix_type, & data_type=new_type) ELSE CALL dbcsr_create(matrix_b, & data_type=new_type, & matrix_type=new_matrix_type, & template=matrix_a) END IF CALL ensure_array_size(matrix_b%index, ub=SIZE(matrix_a%index), & memory_type=dbcsr_get_index_memory_type(matrix_b)) ! ! copy index and data matrix_b%index(1:SIZE(matrix_a%index)) = matrix_a%index(:) IF (.NOT. shallow) THEN IF (matrix_a%nze > dbcsr_get_data_size(matrix_a)) & DBCSR_ABORT("Source matrix sizes not consistent!") CALL dbcsr_data_ensure_size(matrix_b%data_area, & dbcsr_data_get_size_referenced(matrix_a%data_area)) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_get_data_type(matrix_b)) & THEN CALL dbcsr_data_copyall(matrix_b%data_area, & matrix_a%data_area) ELSE CALL dbcsr_data_convert(matrix_b%data_area, & matrix_a%data_area, drop_real=keep_imaginary) END IF ELSE IF (dbcsr_get_data_type(matrix_a) .NE. dbcsr_get_data_type(matrix_b)) & DBCSR_ABORT("Shallow copy only possible when retaining data type.") CALL dbcsr_switch_data_area(matrix_b, matrix_a%data_area) END IF ! ! the row_p, col_i and blk_p ... CALL dbcsr_repoint_index(matrix_b) matrix_b%nze = matrix_a%nze matrix_b%nblks = matrix_b%nblks matrix_b%valid = .TRUE. matrix_b%sparsity_id = matrix_a%sparsity_id END IF CALL timestop(handle) END SUBROUTINE dbcsr_copy SUBROUTINE dbcsr_copy_into_existing(matrix_b, matrix_a) !! copy a matrix, retaining current sparsity TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b !! target DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! source DBCSR matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_into_existing' INTEGER :: col_size, data_type, dst_col, dst_row, & handle, rel, row_size, src_col, & src_cs, src_row, src_rs LOGICAL :: dst_tr, making_symmetric, neg_imag, & neg_real, src_tr TYPE(dbcsr_data_obj) :: dst_data, src_data TYPE(dbcsr_iterator) :: dst_iter, src_iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (.NOT. dbcsr_valid_index(matrix_b)) & DBCSR_ABORT("Matrix_b is not valid") IF (dbcsr_get_data_type(matrix_b) .NE. dbcsr_get_data_type(matrix_a)) & DBCSR_ABORT("Matrices have different data types.") data_type = dbcsr_get_data_type(matrix_b) neg_real = matrix_b%negate_real neg_imag = matrix_b%negate_imaginary making_symmetric = dbcsr_has_symmetry(matrix_b) & .AND. .NOT. dbcsr_has_symmetry(matrix_a) IF (making_symmetric) THEN CALL dbcsr_copy_into_existing_sym(matrix_b, matrix_a) CALL timestop(handle) RETURN END IF CALL dbcsr_data_init(src_data) CALL dbcsr_data_init(dst_data) CALL dbcsr_data_new(src_data, data_type) CALL dbcsr_data_new(dst_data, data_type) CALL dbcsr_iterator_start(src_iter, matrix_a) CALL dbcsr_iterator_start(dst_iter, matrix_b) ! Iterate through the blocks of the source and destination ! matrix. There are three possibilities: 1. copy the data for ! blocks present in both; 2 skip source blocks not present in the ! target; 3 zero blocks not present in the source. IF (dbcsr_iterator_blocks_left(src_iter)) THEN CALL dbcsr_iterator_next_block(src_iter, src_row, src_col, src_data, & src_tr) ELSE src_row = 0; src_col = 0 END IF DO WHILE (dbcsr_iterator_blocks_left(dst_iter)) CALL dbcsr_iterator_next_block(dst_iter, dst_row, dst_col, dst_data, & dst_tr, row_size=row_size, col_size=col_size) ! Now find the source position that is greater or equal to the ! target one. I.e, skip blocks that the target doesn't have. rel = pos_relation(dst_row, dst_col, src_row, src_col) DO WHILE (rel .EQ. 1 .AND. dbcsr_iterator_blocks_left(src_iter)) CALL dbcsr_iterator_next_block(src_iter, src_row, src_col, & src_data, src_tr, row_size=src_rs, col_size=src_cs) rel = pos_relation(dst_row, dst_col, src_row, src_col) END DO SELECT CASE (rel) CASE (-1, 1) ! Target lags source or ran out of source CALL dbcsr_data_clear(dst_data) CASE (0) ! Copy the data IF (dbcsr_data_get_size(src_data) .NE. dbcsr_data_get_size(dst_data)) & DBCSR_ABORT("Block sizes not equal!") IF (src_tr .EQV. dst_tr) THEN CALL dbcsr_data_copyall(dst_data, src_data) ELSE CALL dbcsr_block_partial_copy(dst=dst_data, dst_tr=dst_tr, & dst_rs=row_size, dst_cs=col_size, & dst_r_lb=1, dst_c_lb=1, & src=src_data, src_tr=src_tr, & src_rs=src_rs, src_cs=src_cs, & src_r_lb=1, src_c_lb=1, & nrow=row_size, ncol=col_size) IF (neg_real) THEN CALL dbcsr_block_real_neg(dst_data, row_size, col_size) END IF IF (neg_imag) THEN CALL dbcsr_block_conjg(dst_data, row_size, col_size) END IF END IF CASE default DBCSR_ABORT("Trouble syncing iterators") END SELECT END DO CALL dbcsr_iterator_stop(src_iter) CALL dbcsr_iterator_stop(dst_iter) CALL dbcsr_data_clear_pointer(src_data) CALL dbcsr_data_clear_pointer(dst_data) CALL dbcsr_data_release(src_data) CALL dbcsr_data_release(dst_data) CALL timestop(handle) END SUBROUTINE dbcsr_copy_into_existing SUBROUTINE dbcsr_copy_into_existing_sym(matrix_b, matrix_a) !! copy a matrix, retaining current sparsity TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b !! target DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! source DBCSR matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_into_existing_sym' INTEGER :: col_size, data_type, dst_col, dst_row, & handle, row_size, src_col, src_cs, & src_row, src_rs LOGICAL :: dst_tr, found, neg_imag, neg_real, src_tr TYPE(dbcsr_data_obj) :: dst_data, src_data TYPE(dbcsr_iterator) :: dst_iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_b) .NE. dbcsr_get_data_type(matrix_a)) & DBCSR_ABORT("Matrices have different data types.") data_type = dbcsr_get_data_type(matrix_b) IF (.NOT. dbcsr_has_symmetry(matrix_b) .OR. dbcsr_has_symmetry(matrix_a)) & DBCSR_ABORT("Must copy from non-symmetric to symmetric matrix.") neg_real = matrix_b%negate_real neg_imag = matrix_b%negate_imaginary CALL dbcsr_data_init(src_data) CALL dbcsr_data_init(dst_data) CALL dbcsr_data_new(src_data, data_type) CALL dbcsr_data_new(dst_data, data_type) CALL dbcsr_iterator_start(dst_iter, matrix_b) ! Iterate through the blocks of the destination matrix. For each ! one, try to find an appropriate source matrix block and copy it ! into the destination matrix. DO WHILE (dbcsr_iterator_blocks_left(dst_iter)) CALL dbcsr_iterator_next_block(dst_iter, dst_row, dst_col, dst_data, & dst_tr, row_size=row_size, col_size=col_size) src_row = dst_row src_col = dst_col IF (checker_tr(dst_row, dst_col)) & CALL swap(src_row, src_col) CALL dbcsr_get_block_p(matrix_a, src_row, src_col, src_data, src_tr, & found=found, row_size=src_rs, col_size=src_cs) IF (.NOT. found) THEN CALL dbcsr_data_clear(dst_data) ELSE IF (checker_tr(dst_row, dst_col)) THEN src_tr = .NOT. src_tr CALL swap(src_rs, src_cs) END IF CALL dbcsr_block_partial_copy(dst=dst_data, dst_tr=dst_tr, & dst_rs=row_size, dst_cs=col_size, & dst_r_lb=1, dst_c_lb=1, & src=src_data, src_tr=src_tr, & src_rs=src_rs, src_cs=src_cs, & src_r_lb=1, src_c_lb=1, & nrow=row_size, ncol=col_size) IF (neg_real .AND. checker_tr(dst_row, dst_col)) THEN CALL dbcsr_block_real_neg(dst_data, row_size, col_size) END IF IF (neg_imag .AND. checker_tr(dst_row, dst_col)) THEN CALL dbcsr_block_conjg(dst_data, row_size, col_size) END IF END IF END DO CALL dbcsr_iterator_stop(dst_iter) CALL dbcsr_data_clear_pointer(src_data) CALL dbcsr_data_clear_pointer(dst_data) CALL dbcsr_data_release(src_data) CALL dbcsr_data_release(dst_data) CALL timestop(handle) END SUBROUTINE dbcsr_copy_into_existing_sym ELEMENTAL FUNCTION pos_relation(row1, col1, row2, col2) RESULT(relation) !! Determines the relation between two matrix positions. INTEGER, INTENT(IN) :: row1, col1, row2, col2 INTEGER :: relation !! Relation between positions 1 and 2. 0: same -1: pos1 < pos2 1: pos1 > pos2 IF (row1 .LT. row2) THEN relation = -1 ELSEIF (row1 .GT. row2) THEN relation = 1 ELSE ! rows are equal, check column IF (col1 .LT. col2) THEN relation = -1 ELSEIF (col1 .GT. col2) THEN relation = 1 ELSE relation = 0 END IF END IF END FUNCTION pos_relation SUBROUTINE dbcsr_copy_submatrix(matrix_b, matrix_a, name, & block_row_bounds, block_column_bounds, & shallow_data) !! Copy a submatrix. TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b !! target DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! source DBCSR matrix CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name !! name of the new matrix INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: block_row_bounds, block_column_bounds !! rows to extract (array of size 2 holding the lower and upper inclusive bounds) !! columns to extract (array of size 2 holding the lower and upper inclusive bounds) LOGICAL, INTENT(IN), OPTIONAL :: shallow_data !! shallow data copy CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy_submatrix' INTEGER :: blk_p, col, handle, nblocks, new_blk, & old_blk, row INTEGER, ALLOCATABLE, DIMENSION(:) :: blkp_list, col_list, row_list LOGICAL :: shallow, tr TYPE(dbcsr_data_obj) :: data_block TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) IF (PRESENT(shallow_data)) THEN shallow = shallow_data ELSE shallow = .FALSE. END IF ! Verify assumptions. IF (PRESENT(block_row_bounds)) THEN IF (SIZE(block_row_bounds) /= 2) & DBCSR_ABORT("Size of bounds specifier must be 2") END IF IF (PRESENT(block_column_bounds)) THEN IF (SIZE(block_column_bounds) /= 2) & DBCSR_ABORT("Size of bounds specifier must be 2") END IF ! Setup target matrix CALL dbcsr_create(matrix_b, name=name, template=matrix_a) CALL dbcsr_finalize(matrix_b) IF (.NOT. shallow) THEN ! Non-shallow copy uses the standard iterator on the source and ! block put on the target. ! !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (data_block, iter, row, col, tr) & !$OMP SHARED (matrix_a, matrix_b,& !$OMP block_row_bounds, block_column_bounds) CALL dbcsr_work_create(matrix_b, work_mutable=.FALSE.) CALL dbcsr_data_init(data_block) CALL dbcsr_data_new(data_block, dbcsr_get_data_type(matrix_a)) CALL dbcsr_iterator_start(iter, matrix_a, dynamic=.TRUE., & dynamic_byrows=.TRUE.) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr) ! Only keep the block if they are within the specified bounds. IF (PRESENT(block_row_bounds)) THEN IF (row .LT. block_row_bounds(1)) CYCLE IF (row .GT. block_row_bounds(2)) CYCLE END IF IF (PRESENT(block_column_bounds)) THEN IF (col .LT. block_column_bounds(1)) CYCLE IF (col .GT. block_column_bounds(2)) CYCLE END IF CALL dbcsr_put_block(matrix_b, row, col, data_block, transposed=tr) END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_block) CALL dbcsr_data_release(data_block) CALL dbcsr_finalize(matrix_b) !$OMP END PARALLEL ELSE ! For the shallow copy the source matrix data is referenced. CALL dbcsr_switch_data_area(matrix_b, matrix_a%data_area) nblocks = dbcsr_get_num_blocks(matrix_a) ! High estimate. ! Shallow copy goes through source's data blocks and inserts ! the only the ones corresponding to the submatrix specifier ! into the target. Block pointers must remain the same as in ! the source. ALLOCATE (row_list(nblocks), col_list(nblocks), blkp_list(nblocks)) ! CALL dbcsr_iterator_start(iter, matrix_a) new_blk = 1 DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, & blk=old_blk, blk_p=blk_p) ! Only keep the block if they are within the specified bounds. IF (PRESENT(block_row_bounds)) THEN IF (row .LT. block_row_bounds(1)) CYCLE IF (row .GT. block_row_bounds(2)) CYCLE END IF IF (PRESENT(block_column_bounds)) THEN IF (col .LT. block_column_bounds(1)) CYCLE IF (col .GT. block_column_bounds(2)) CYCLE END IF row_list(new_blk) = row col_list(new_blk) = col blkp_list(new_blk) = blk_p new_blk = new_blk + 1 END DO new_blk = new_blk - 1 CALL dbcsr_iterator_stop(iter) CALL dbcsr_reserve_blocks(matrix_b, row_list(1:new_blk), & col_list(1:new_blk), blkp_list(1:new_blk)) END IF ! CALL timestop(handle) END SUBROUTINE dbcsr_copy_submatrix SUBROUTINE dbcsr_crop_matrix(matrix_b, matrix_a, & full_row_bounds, full_column_bounds, & shallow_data) !! Crop and copies a matrix. TYPE(dbcsr_type), INTENT(INOUT) :: matrix_b !! target DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! source DBCSR matrix INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: full_row_bounds, full_column_bounds !! rows to extract (array of size 2 holding the lower and upper inclusive bounds) !! columns to extract (array of size 2 holding the lower and upper inclusive bounds) LOGICAL, INTENT(IN), OPTIONAL :: shallow_data CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_crop_matrix' INTEGER :: col, f_col_f, f_row_f, handle, l_col_l, & l_row_l, row INTEGER, DIMENSION(2) :: block_col_bounds, block_row_bounds LOGICAL :: part_col, part_f_col, part_f_row, & part_l_col, part_l_row, part_row, & shallow, tr TYPE(dbcsr_data_obj) :: data_block TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) part_l_col = .FALSE. part_f_col = .FALSE. part_l_row = .FALSE. part_f_row = .FALSE. IF (PRESENT(shallow_data)) THEN shallow = shallow_data ELSE shallow = .FALSE. END IF block_row_bounds = 0 block_col_bounds = 0 part_col = .FALSE. part_row = .FALSE. ! ! If row bounds are present, they must be converted to block ! addressing. IF (PRESENT(full_row_bounds)) THEN IF (SIZE(full_row_bounds) /= 2) & DBCSR_ABORT("Size of bounds specifier must be 2") IF (full_row_bounds(1) < 0) & DBCSR_ABORT("Invalid first row bound.") IF (full_row_bounds(2) > dbcsr_nfullrows_total(matrix_a)) & DBCSR_ABORT("Invalid last row bound.") IF (full_row_bounds(1) .EQ. 0) THEN block_row_bounds(1) = 1 ELSE CALL find_block_of_element(full_row_bounds(1), block_row_bounds(1), & dbcsr_nblkrows_total(matrix_a), & dbcsr_row_block_offsets(matrix_a), & hint=0) part_f_row = array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(1)) & .NE. full_row_bounds(1) END IF f_row_f = -7 IF (part_f_row) THEN ! Block offset of last cleared row f_row_f = full_row_bounds(1) - & array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(1)) END IF IF (full_row_bounds(2) .EQ. 0) THEN block_row_bounds(2) = dbcsr_nblkrows_total(matrix_a) ELSE CALL find_block_of_element(full_row_bounds(2), block_row_bounds(2), & dbcsr_nblkrows_total(matrix_a), & dbcsr_row_block_offsets(matrix_a), & hint=0) part_l_row = array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(2) + 1) - 1 & .NE. full_row_bounds(2) END IF ! Block offset of first cleared row l_row_l = -7 IF (part_l_row) THEN l_row_l = 2 + full_row_bounds(2) - & array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(2)) END IF part_row = part_f_row .OR. part_l_row END IF ! ! If column bounds are present, they must be converted to block ! addressing. IF (PRESENT(full_column_bounds)) THEN IF (SIZE(full_column_bounds) /= 2) & DBCSR_ABORT("Size of bounds specifier must be 2") IF (full_column_bounds(1) < 0) & DBCSR_ABORT("Invalid first column bound.") IF (full_column_bounds(2) > dbcsr_nfullcols_total(matrix_a)) & DBCSR_ABORT("Invalid last column bound.") IF (full_column_bounds(1) .EQ. 0) THEN block_col_bounds(1) = 1 ELSE CALL find_block_of_element(full_column_bounds(1), block_col_bounds(1), & dbcsr_nblkcols_total(matrix_a), & dbcsr_col_block_offsets(matrix_a), & hint=0) part_f_col = array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(1)) & .NE. full_column_bounds(1) END IF f_col_f = -7 IF (part_f_col) THEN ! Block offset of last cleared column f_col_f = full_column_bounds(1) - & array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(1)) END IF IF (full_column_bounds(2) .EQ. 0) THEN block_col_bounds(2) = dbcsr_nblkcols_total(matrix_a) ELSE CALL find_block_of_element(full_column_bounds(2), block_col_bounds(2), & dbcsr_nblkcols_total(matrix_a), & dbcsr_col_block_offsets(matrix_a), & hint=0) part_l_col = array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(2) + 1) - 1 & .NE. full_column_bounds(2) END IF l_col_l = -7 IF (part_l_col) THEN ! Block offset of first cleared column l_col_l = 2 + full_column_bounds(2) - & array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(2)) END IF part_col = part_f_col .OR. part_l_col END IF ! ! First copy the blocks then perform the intra-block zeroing. CALL dbcsr_copy_submatrix(matrix_b, matrix_a, & block_row_bounds=block_row_bounds, & block_column_bounds=block_col_bounds, & shallow_data=shallow) IF (part_row .OR. part_col) THEN !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (data_block, iter, row, col, tr) & !$OMP SHARED (matrix_b,& !$OMP part_row, part_f_row, part_l_row, f_row_f, l_row_l, & !$OMP part_col, part_f_col, part_l_col, f_col_f, l_col_l,& !$OMP block_row_bounds, block_col_bounds) CALL dbcsr_data_init(data_block) CALL dbcsr_data_new(data_block, dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix_b))) CALL dbcsr_iterator_start(iter, matrix_b, & dynamic=.TRUE., dynamic_byrows=.TRUE.) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr) IF (part_row) THEN IF (row .LT. block_row_bounds(1)) CYCLE IF (row .GT. block_row_bounds(2)) CYCLE END IF IF (part_col) THEN IF (col .LT. block_col_bounds(1)) CYCLE IF (col .GT. block_col_bounds(2)) CYCLE END IF IF (part_row) THEN IF (part_f_row .AND. row .EQ. block_row_bounds(1)) THEN CALL dbcsr_data_clear(data_block, ub=f_row_f, tr=tr) END IF IF (part_l_row .AND. row .EQ. block_row_bounds(2)) THEN CALL dbcsr_data_clear(data_block, lb=l_row_l, tr=tr) END IF END IF IF (part_col) THEN IF (part_f_col .AND. col .EQ. block_col_bounds(1)) THEN CALL dbcsr_data_clear(data_block, ub2=f_col_f, tr=tr) END IF IF (part_l_col .AND. col .EQ. block_col_bounds(2)) THEN CALL dbcsr_data_clear(data_block, lb2=l_col_l, tr=tr) END IF END IF END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_block) CALL dbcsr_data_release(data_block) CALL dbcsr_finalize(matrix_b) !$OMP END PARALLEL END IF ! CALL timestop(handle) END SUBROUTINE dbcsr_crop_matrix SUBROUTINE dbcsr_triu(matrix_a) !! triu of a dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! the matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_triu' INTEGER :: blk, blk_nze, col, col_size, handle, i, & j, row, row_size LOGICAL :: tr REAL(dp), DIMENSION(:, :), POINTER :: DATA TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) CALL dbcsr_iterator_start(iter, matrix_a) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, DATA, tr, & block_number=blk, row_size=row_size, col_size=col_size) blk_nze = row_size*col_size IF (row .GT. col) CALL dbcsr_remove_block(matrix_a, row, col, blk_nze, blk) IF (row .EQ. col) THEN DO j = 1, col_size DO i = j + 1, row_size DATA(i, j) = 0.0_dp END DO END DO END IF END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_finalize(matrix_a) CALL timestop(handle) END SUBROUTINE dbcsr_triu SUBROUTINE dbcsr_filter_anytype(matrix, eps, method, & use_absolute, filter_diag) !! filter a dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! the matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: eps !! the threshold INTEGER, INTENT(IN), OPTIONAL :: method !! how the matrix is filtered LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag !! NYI CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_filter_anytype' COMPLEX(KIND=real_4), DIMENSION(:), POINTER :: data_c COMPLEX(KIND=real_8), DIMENSION(:), POINTER :: data_z INTEGER :: blk, blk_nze, col, col_size, handle, & my_method, row, row_size, data_type LOGICAL :: gt0, my_filter_diag, tr REAL(KIND=real_4) :: nrm_s REAL(KIND=real_4), DIMENSION(:), POINTER :: data_s REAL(KIND=real_8) :: my_absolute, nrm_d REAL(KIND=real_8), DIMENSION(:), POINTER :: data_d TYPE(dbcsr_iterator) :: iter REAL(KIND=real_8), EXTERNAL :: DZNRM2 #if defined (__ACCELERATE) REAL(KIND=real_8), EXTERNAL :: SCNRM2 #else REAL(KIND=real_4), EXTERNAL :: SCNRM2 #endif ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) my_method = dbcsr_filter_frobenius IF (PRESENT(method)) my_method = method my_absolute = 1.0_dp IF (PRESENT(use_absolute)) my_absolute = dbcsr_maxabs(matrix) my_filter_diag = .TRUE. IF (PRESENT(filter_diag)) my_filter_diag = filter_diag SELECT CASE (eps%data_type) CASE (dbcsr_type_real_4) gt0 = eps%r_sp .GT. 0.0_real_4 CASE (dbcsr_type_real_8) gt0 = eps%r_dp .GT. 0.0_real_8 CASE (dbcsr_type_complex_4) gt0 = ABS(eps%c_sp) .GT. 0.0_real_4 CASE (dbcsr_type_complex_8) gt0 = ABS(eps%c_dp) .GT. 0.0_real_8 CASE default gt0 = .FALSE. END SELECT IF (gt0) THEN data_type = dbcsr_get_data_type(matrix) !$OMP PARALLEL DEFAULT(NONE) PRIVATE(iter,row,col,data_s,data_d,data_c,data_z,tr, & !$OMP blk,row_size,col_size,blk_nze,nrm_d,nrm_s) & !$OMP SHARED(my_method,my_absolute,eps,matrix,data_type) CALL dbcsr_iterator_start(iter, matrix, contiguous_pointers=.TRUE.) DO WHILE (dbcsr_iterator_blocks_left(iter)) SELECT CASE (data_type) CASE (dbcsr_type_real_4) CALL dbcsr_iterator_next_block(iter, row, col, data_s, tr, blk, & row_size, col_size) blk_nze = row_size*col_size IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks SELECT CASE (my_method) CASE (dbcsr_filter_frobenius) ! ! Frobenius based nrm_s = norm2(data_s) IF (nrm_s .LT. my_absolute*eps%r_sp) & CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) CASE DEFAULT DBCSR_ABORT("Only Frobenius based filtering") END SELECT CASE (dbcsr_type_real_8) CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk, & row_size, col_size) blk_nze = row_size*col_size IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks SELECT CASE (my_method) CASE (dbcsr_filter_frobenius) ! ! Frobenius based nrm_d = norm2(data_d) IF (nrm_d .LT. my_absolute*eps%r_dp) & CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) CASE DEFAULT DBCSR_ABORT("Only Frobenius based filtering") END SELECT CASE (dbcsr_type_complex_4) CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk, & row_size, col_size) blk_nze = row_size*col_size IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks SELECT CASE (my_method) CASE (dbcsr_filter_frobenius) ! ! Frobenius based nrm_d = SCNRM2(SIZE(data_c), data_c(1), 1) IF (nrm_d .LT. my_absolute*eps%r_dp) & CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) CASE DEFAULT DBCSR_ABORT("Only Frobenius based filtering") END SELECT CASE (dbcsr_type_complex_8) CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk, & row_size, col_size) blk_nze = row_size*col_size IF (blk_nze .EQ. 0) CYCLE ! Skip empty blocks SELECT CASE (my_method) CASE (dbcsr_filter_frobenius) ! ! Frobenius based nrm_d = DZNRM2(SIZE(data_z), data_z(1), 1) IF (nrm_d .LT. my_absolute*eps%r_dp) & CALL dbcsr_remove_block(matrix, row, col, blk_nze, blk) CASE DEFAULT DBCSR_ABORT("Only Frobenius based filtering") END SELECT CASE DEFAULT DBCSR_ABORT("Wrong data type") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_finalize(matrix, reshuffle=.TRUE.) !$OMP END PARALLEL CALL dbcsr_index_compact(matrix) END IF CALL timestop(handle) END SUBROUTINE dbcsr_filter_anytype SUBROUTINE dbcsr_norm_scalar(matrix, which_norm, norm_scalar) !! compute a norm of a dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! the matrix INTEGER, INTENT(IN) :: which_norm REAL(KIND=real_8), INTENT(OUT) :: norm_scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_norm_scalar' INTEGER :: handle ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) SELECT CASE (which_norm) CASE (dbcsr_norm_frobenius) norm_scalar = dbcsr_frobenius_norm(matrix) CASE (dbcsr_norm_maxabsnorm) norm_scalar = dbcsr_maxabs(matrix) CASE (dbcsr_norm_gershgorin) norm_scalar = dbcsr_gershgorin_norm(matrix) CASE DEFAULT DBCSR_ABORT("this norm is NYI") END SELECT CALL timestop(handle) END SUBROUTINE SUBROUTINE dbcsr_norm_r8_vec(matrix, which_norm, norm_vector) TYPE(dbcsr_type), INTENT(INOUT) :: matrix INTEGER, INTENT(IN) :: which_norm REAL(KIND=real_8), DIMENSION(:), INTENT(OUT), & TARGET, CONTIGUOUS :: norm_vector REAL(KIND=real_8), DIMENSION(:), POINTER, CONTIGUOUS :: v_p TYPE(dbcsr_data_obj) :: norm_vector_a CALL dbcsr_data_init(norm_vector_a) CALL dbcsr_data_new(norm_vector_a, dbcsr_type_real_8) v_p => norm_vector CALL dbcsr_data_set_pointer(norm_vector_a, v_p) CALL dbcsr_norm_vec(matrix, which_norm, norm_vector_a) CALL dbcsr_data_clear_pointer(norm_vector_a) CALL dbcsr_data_release(norm_vector_a) END SUBROUTINE dbcsr_norm_r8_vec SUBROUTINE dbcsr_norm_vec(matrix, which_norm, norm_vector) !! compute the column norms of the dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! the matrix INTEGER, INTENT(IN) :: which_norm TYPE(dbcsr_data_obj), INTENT(INOUT) :: norm_vector CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_norm_vec' INTEGER :: blk, col, col_offset, i, j, row, & row_offset, handle LOGICAL :: tr TYPE(dbcsr_data_obj) :: data_a TYPE(dbcsr_iterator) :: iter CALL timeset(routineN, handle) SELECT CASE (which_norm) CASE (dbcsr_norm_column) IF (dbcsr_data_get_type(norm_vector) /= dbcsr_get_data_type(matrix)) & DBCSR_ABORT("Mismatched vector/matrix data types") IF (dbcsr_has_symmetry(matrix)) THEN IF (dbcsr_data_get_size(norm_vector) < dbcsr_nfullrows_total(matrix)) & DBCSR_ABORT("Passed vector too small") END IF IF (dbcsr_data_get_size(norm_vector) < dbcsr_nfullcols_total(matrix)) & DBCSR_ABORT("Passed vector too small") CALL dbcsr_data_init(data_a) CALL dbcsr_data_new(data_a, dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix))) CALL dbcsr_data_clear(norm_vector) CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, data_a, tr, & blk, row_offset=row_offset, col_offset=col_offset) SELECT CASE (dbcsr_get_data_type(matrix)) CASE (dbcsr_type_real_4) IF (dbcsr_has_symmetry(matrix) .AND. row .NE. col) THEN DO j = 1, SIZE(data_a%d%r2_sp, 2) DO i = 1, SIZE(data_a%d%r2_sp, 1) norm_vector%d%r_sp(col_offset + j - 1) & = norm_vector%d%r_sp(col_offset + j - 1) & + data_a%d%r2_sp(i, j)**2 norm_vector%d%r_sp(row_offset + i - 1) & = norm_vector%d%r_sp(row_offset + i - 1) & + data_a%d%r2_sp(i, j)**2 END DO END DO ELSE DO j = 1, SIZE(data_a%d%r2_sp, 2) DO i = 1, SIZE(data_a%d%r2_sp, 1) norm_vector%d%r_sp(col_offset + j - 1) & = norm_vector%d%r_sp(col_offset + j - 1) & + data_a%d%r2_sp(i, j)*data_a%d%r2_sp(i, j) END DO END DO END IF CASE (dbcsr_type_real_8) IF (dbcsr_has_symmetry(matrix) .AND. row .NE. col) THEN DO j = 1, SIZE(data_a%d%r2_dp, 2) DO i = 1, SIZE(data_a%d%r2_dp, 1) norm_vector%d%r_dp(col_offset + j - 1) & = norm_vector%d%r_dp(col_offset + j - 1) & + data_a%d%r2_dp(i, j)**2 norm_vector%d%r_dp(row_offset + i - 1) & = norm_vector%d%r_dp(row_offset + i - 1) & + data_a%d%r2_dp(i, j)**2 END DO END DO ELSE DO j = 1, SIZE(data_a%d%r2_dp, 2) DO i = 1, SIZE(data_a%d%r2_dp, 1) norm_vector%d%r_dp(col_offset + j - 1) & = norm_vector%d%r_dp(col_offset + j - 1) & + data_a%d%r2_dp(i, j)*data_a%d%r2_dp(i, j) END DO END DO END IF CASE DEFAULT DBCSR_ABORT("Only real values") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL dbcsr_data_clear_pointer(data_a) CALL dbcsr_data_release(data_a) SELECT CASE (dbcsr_get_data_type(matrix)) CASE (dbcsr_type_real_4) CALL mp_sum(norm_vector%d%r_sp, & dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) norm_vector%d%r_sp = SQRT(norm_vector%d%r_sp) CASE (dbcsr_type_real_8) CALL mp_sum(norm_vector%d%r_dp, & dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) norm_vector%d%r_dp = SQRT(norm_vector%d%r_dp) END SELECT CASE DEFAULT DBCSR_ABORT("this norm is NYI") END SELECT CALL timestop(handle) END SUBROUTINE dbcsr_norm_vec FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm) !! compute a norm of a dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! the matrix REAL(KIND=real_8) :: norm CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_gershgorin_norm' COMPLEX(KIND=real_4), DIMENSION(:, :), POINTER :: data_c COMPLEX(KIND=real_8), DIMENSION(:, :), POINTER :: data_z INTEGER :: blk, col, col_offset, handle, i, j, nc, & nr, row, row_offset LOGICAL :: any_sym, tr REAL(KIND=real_4), DIMENSION(:, :), POINTER :: data_r REAL(KIND=real_8), DIMENSION(:, :), POINTER :: data_d REAL(real_8), ALLOCATABLE, DIMENSION(:) :: buff_d TYPE(dbcsr_iterator) :: iter CALL timeset(routineN, handle) nr = dbcsr_nfullrows_total(matrix) nc = dbcsr_nfullcols_total(matrix) any_sym = dbcsr_get_matrix_type(matrix) .EQ. dbcsr_type_symmetric .OR. & dbcsr_get_matrix_type(matrix) .EQ. dbcsr_type_antisymmetric IF (nr .NE. nc) & DBCSR_ABORT("not a square matrix") norm = 0.0_dp ALLOCATE (buff_d(nr)) buff_d = 0.0_dp CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) SELECT CASE (dbcsr_get_data_type(matrix)) CASE (dbcsr_type_real_4) CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk, & row_offset=row_offset, col_offset=col_offset) DO j = 1, SIZE(data_r, 2) DO i = 1, SIZE(data_r, 1) buff_d(row_offset + i - 1) = buff_d(row_offset + i - 1) + ABS(data_r(i, j)) IF (any_sym .AND. row .NE. col) & buff_d(col_offset + j - 1) = buff_d(col_offset + j - 1) + ABS(data_r(i, j)) END DO END DO CASE (dbcsr_type_real_8) CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk, & row_offset=row_offset, col_offset=col_offset) DO j = 1, SIZE(data_d, 2) DO i = 1, SIZE(data_d, 1) buff_d(row_offset + i - 1) = buff_d(row_offset + i - 1) + ABS(data_d(i, j)) IF (any_sym .AND. row .NE. col) & buff_d(col_offset + j - 1) = buff_d(col_offset + j - 1) + ABS(data_d(i, j)) END DO END DO CASE (dbcsr_type_complex_4) CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk, & row_offset=row_offset, col_offset=col_offset) DO j = 1, SIZE(data_c, 2) DO i = 1, SIZE(data_c, 1) buff_d(row_offset + i - 1) = buff_d(row_offset + i - 1) + ABS(data_c(i, j)) IF (any_sym .AND. row .NE. col) & DBCSR_ABORT("Only nonsymmetric matrix so far") ! buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_c(i,j)) END DO END DO CASE (dbcsr_type_complex_8) CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk, & row_offset=row_offset, col_offset=col_offset) DO j = 1, SIZE(data_z, 2) DO i = 1, SIZE(data_z, 1) buff_d(row_offset + i - 1) = buff_d(row_offset + i - 1) + ABS(data_z(i, j)) IF (any_sym .AND. row .NE. col) & DBCSR_ABORT("Only nonsymmetric matrix so far") ! buff_d(col_offset+j-1) = buff_d(col_offset+j-1) + ABS(data_z(i,j)) END DO END DO CASE DEFAULT DBCSR_ABORT("Wrong data type") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL mp_sum(buff_d, dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) norm = MAXVAL(buff_d) DEALLOCATE (buff_d) CALL timestop(handle) END FUNCTION dbcsr_gershgorin_norm FUNCTION dbcsr_maxabs(matrix) RESULT(norm) !! compute a norm of a dbcsr matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! the matrix REAL(real_8) :: norm COMPLEX(KIND=real_4), DIMENSION(:, :), POINTER :: data_c COMPLEX(KIND=real_8), DIMENSION(:, :), POINTER :: data_z INTEGER :: blk, col, row LOGICAL :: tr REAL(KIND=real_4), DIMENSION(:, :), POINTER :: data_r REAL(KIND=real_8), DIMENSION(:, :), POINTER :: data_d TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- norm = 0.0_dp CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) SELECT CASE (dbcsr_get_data_type(matrix)) CASE (dbcsr_type_real_4) CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk) norm = MAX(norm, REAL(MAXVAL(ABS(data_r)), dp)) CASE (dbcsr_type_real_8) CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk) norm = MAX(norm, MAXVAL(ABS(data_d))) CASE (dbcsr_type_complex_4) CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk) norm = MAX(norm, REAL(MAXVAL(ABS(data_c)), dp)) CASE (dbcsr_type_complex_8) CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk) norm = MAX(norm, MAXVAL(ABS(data_z))) CASE DEFAULT DBCSR_ABORT("Wrong data type") END SELECT END DO CALL dbcsr_iterator_stop(iter) CALL mp_max(norm, dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) END FUNCTION dbcsr_maxabs FUNCTION dbcsr_frobenius_norm(matrix, local) RESULT(norm) !! compute a norm of a dbcsr matrix TYPE(dbcsr_type), INTENT(IN) :: matrix !! the matrix LOGICAL, INTENT(in), OPTIONAL :: local REAL(KIND=real_8) :: norm CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_frobenius_norm' COMPLEX(KIND=real_4), DIMENSION(:, :), POINTER :: data_c COMPLEX(KIND=real_8), DIMENSION(:, :), POINTER :: data_z INTEGER :: blk, col, handle, row LOGICAL :: any_sym, my_local, tr REAL(KIND=real_4), DIMENSION(:, :), POINTER :: data_r REAL(KIND=real_8), DIMENSION(:, :), POINTER :: data_d REAL(real_8) :: fac TYPE(dbcsr_iterator) :: iter ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) my_local = .FALSE. IF (PRESENT(local)) my_local = local any_sym = dbcsr_get_matrix_type(matrix) .EQ. dbcsr_type_symmetric .OR. & dbcsr_get_matrix_type(matrix) .EQ. dbcsr_type_antisymmetric norm = 0.0_dp CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) SELECT CASE (dbcsr_get_data_type(matrix)) CASE (dbcsr_type_real_4) CALL dbcsr_iterator_next_block(iter, row, col, data_r, tr, blk) fac = 1.0_dp IF (any_sym .AND. row .NE. col) fac = 2.0_dp norm = norm + fac*SUM(data_r**2) CASE (dbcsr_type_real_8) CALL dbcsr_iterator_next_block(iter, row, col, data_d, tr, blk) fac = 1.0_dp IF (any_sym .AND. row .NE. col) fac = 2.0_dp norm = norm + fac*SUM(data_d**2) CASE (dbcsr_type_complex_4) CALL dbcsr_iterator_next_block(iter, row, col, data_c, tr, blk) fac = 1.0_dp IF (any_sym .AND. row .NE. col) & DBCSR_ABORT("Only nonsymmetric matrix so far") norm = norm + fac*REAL(SUM(CONJG(data_c)*data_c), KIND=real_8) CASE (dbcsr_type_complex_8) CALL dbcsr_iterator_next_block(iter, row, col, data_z, tr, blk) fac = 1.0_dp IF (any_sym .AND. row .NE. col) & DBCSR_ABORT("Only nonsymmetric matrix so far") norm = norm + fac*REAL(SUM(CONJG(data_z)*data_z), KIND=real_8) CASE DEFAULT DBCSR_ABORT("Wrong data type") END SELECT END DO CALL dbcsr_iterator_stop(iter) IF (.NOT. my_local) CALL mp_sum(norm, dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) norm = SQRT(norm) CALL timestop(handle) END FUNCTION dbcsr_frobenius_norm SUBROUTINE dbcsr_sum_replicated(matrix) !! Sums blocks in a replicated dbcsr matrix, which has the same structure on all ranks. TYPE(dbcsr_type), INTENT(inout) :: matrix !! dbcsr matrix to operate on CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_sum_replicated' INTEGER :: handle, index_checksum, mynode, & numnodes INTEGER, ALLOCATABLE, DIMENSION(:) :: all_checksums TYPE(dbcsr_mp_obj) :: mp TYPE(mp_comm_type) :: comm ! --------------------------------------------------------------------------- CALL timeset(routineN, handle) mp = dbcsr_distribution_mp(dbcsr_distribution(matrix)) comm = dbcsr_mp_group(mp) numnodes = dbcsr_mp_numnodes(mp) mynode = dbcsr_mp_mynode(mp) ! ALLOCATE (all_checksums(numnodes)) CALL dbcsr_index_checksum(matrix, index_checksum) CALL mp_allgather(index_checksum, all_checksums, comm) ! IF (.NOT. ALL(all_checksums .EQ. index_checksum)) & DBCSR_ABORT("Replicated matrices do not all have the same index structure.") ! SELECT CASE (dbcsr_data_get_type(matrix%data_area)) CASE (dbcsr_type_real_4) CALL mp_sum(matrix%data_area%d%r_sp, comm) CASE (dbcsr_type_real_8) CALL mp_sum(matrix%data_area%d%r_dp, comm) CASE (dbcsr_type_complex_4) CALL mp_sum(matrix%data_area%d%c_sp, comm) CASE (dbcsr_type_complex_8) CALL mp_sum(matrix%data_area%d%c_dp, comm) CASE default DBCSR_ABORT("Incorrect data type") END SELECT ! CALL timestop(handle) END SUBROUTINE dbcsr_sum_replicated FUNCTION dbcsr_block_in_limits(row, col, block_row_limits, block_column_limits) !! check if a block is not in the limits INTEGER, INTENT(in) :: row, col INTEGER, DIMENSION(2), INTENT(in), OPTIONAL :: block_row_limits, block_column_limits LOGICAL :: dbcsr_block_in_limits dbcsr_block_in_limits = .TRUE. IF (PRESENT(block_row_limits)) THEN IF (row .LT. block_row_limits(1)) dbcsr_block_in_limits = .FALSE. IF (row .GT. block_row_limits(2)) dbcsr_block_in_limits = .FALSE. END IF IF (PRESENT(block_column_limits)) THEN IF (col .LT. block_column_limits(1)) dbcsr_block_in_limits = .FALSE. IF (col .GT. block_column_limits(2)) dbcsr_block_in_limits = .FALSE. END IF END FUNCTION dbcsr_block_in_limits SUBROUTINE dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, & nfullrows_total, nfullcols_total, & nblkrows_local, nblkcols_local, & nfullrows_local, nfullcols_local, & my_prow, my_pcol, & local_rows, local_cols, proc_row_dist, proc_col_dist, & row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, data_area, & matrix_type, data_type, group) !! Gets information about a matrix TYPE(dbcsr_type), INTENT(IN) :: matrix !! matrix to query INTEGER, INTENT(OUT), OPTIONAL :: nblkrows_total, nblkcols_total, nfullrows_total, & nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, & my_prow, my_pcol INTEGER, DIMENSION(:), OPTIONAL, POINTER :: local_rows, local_cols, proc_row_dist, & proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset TYPE(dbcsr_distribution_obj), INTENT(OUT), & OPTIONAL :: distribution !! the data distribution of the matrix CHARACTER(len=*), INTENT(OUT), OPTIONAL :: name !! matrix name TYPE(dbcsr_data_obj), INTENT(OUT), OPTIONAL :: data_area !! data_area CHARACTER, OPTIONAL :: matrix_type !! matrix type (regular, symmetric, see dbcsr_types.F for values) INTEGER, OPTIONAL :: data_type !! data type (single/double precision real/complex) TYPE(mp_comm_type), INTENT(OUT), OPTIONAL :: group ! --------------------------------------------------------------------------- !vw avoid massive printing of warnings !DBCSR_WARN("Invalid matrix") IF (PRESENT(nblkrows_total)) nblkrows_total = matrix%nblkrows_total IF (PRESENT(nblkcols_total)) nblkcols_total = matrix%nblkcols_total IF (PRESENT(nfullrows_total)) nfullrows_total = matrix%nfullrows_total IF (PRESENT(nfullcols_total)) nfullcols_total = matrix%nfullcols_total IF (PRESENT(nblkrows_local)) nblkrows_local = matrix%nblkrows_local IF (PRESENT(nblkcols_local)) nblkcols_local = matrix%nblkcols_local IF (PRESENT(nfullrows_local)) nfullrows_local = matrix%nfullrows_local IF (PRESENT(nfullcols_local)) nfullcols_local = matrix%nfullcols_local IF (PRESENT(row_blk_size)) row_blk_size => array_data(matrix%row_blk_size) IF (PRESENT(col_blk_size)) col_blk_size => array_data(matrix%col_blk_size) IF (PRESENT(row_blk_offset)) row_blk_offset => array_data(matrix%row_blk_offset) IF (PRESENT(col_blk_offset)) col_blk_offset => array_data(matrix%col_blk_offset) IF (PRESENT(distribution)) distribution = matrix%dist IF (PRESENT(name)) name = matrix%name IF (PRESENT(data_area)) data_area = matrix%data_area IF (PRESENT(data_type)) data_type = matrix%data_type IF (PRESENT(local_rows)) local_rows => dbcsr_distribution_local_rows(matrix%dist) IF (PRESENT(local_cols)) local_cols => dbcsr_distribution_local_cols(matrix%dist) IF (PRESENT(proc_row_dist)) proc_row_dist => dbcsr_distribution_row_dist(matrix%dist) IF (PRESENT(proc_col_dist)) proc_col_dist => dbcsr_distribution_col_dist(matrix%dist) IF (PRESENT(my_prow)) my_prow = dbcsr_mp_myprow(dbcsr_distribution_mp(matrix%dist)) IF (PRESENT(my_pcol)) my_pcol = dbcsr_mp_mypcol(dbcsr_distribution_mp(matrix%dist)) IF (PRESENT(matrix_type)) matrix_type = dbcsr_get_matrix_type(matrix) IF (PRESENT(group)) group = dbcsr_mp_group(matrix%dist%d%mp_env) ! a shortcut !IF (PRESENT(matrix_type)) THEN ! matrix_type = dbcsr_get_matrix_type(matrix) ! IF (matrix_type .EQ. dbcsr_type_invalid) & ! DBCSR_ABORT("Incorrect symmetry") !ENDIF END SUBROUTINE dbcsr_get_info FUNCTION dbcsr_may_be_dense(matrix, occ_thresh) RESULT(may_be_dense) !! Returns whether the matrix could be represented in a dense form TYPE(dbcsr_type), INTENT(IN) :: matrix !! matrix REAL(real_8), INTENT(in) :: occ_thresh LOGICAL :: may_be_dense !! use the mutable and not append-only working structures REAL(real_8) :: occ ! --------------------------------------------------------------------------- occ = dbcsr_get_occupation(matrix) may_be_dense = .NOT. (occ .LT. occ_thresh) ! make sure every proc sees the same CALL mp_sum(may_be_dense, dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) END FUNCTION dbcsr_may_be_dense FUNCTION dbcsr_get_occupation(matrix) RESULT(occupation) !! Returns the occupation of the matrix TYPE(dbcsr_type), INTENT(IN) :: matrix !! matrix from which to get the occupation REAL(KIND=real_8) :: occupation INTEGER :: nfullcols, nfullrows INTEGER(KIND=int_8) :: nze_global INTEGER, DIMENSION(:), POINTER :: row_blk_size nze_global = matrix%nze CALL mp_sum(nze_global, dbcsr_mp_group(dbcsr_distribution_mp(matrix%dist))) nfullrows = dbcsr_nfullrows_total(matrix) nfullcols = dbcsr_nfullcols_total(matrix) row_blk_size => array_data(matrix%row_blk_size) IF (nfullrows .NE. 0 .AND. nfullcols .NE. 0) THEN IF (dbcsr_has_symmetry(matrix)) THEN IF (2*nze_global .EQ. & (INT(nfullrows, KIND=int_8)*INT(nfullrows + 1, KIND=int_8) + SUM(row_blk_size*(row_blk_size - 1)))) THEN occupation = 1.0_real_8 ELSE occupation = 2.0_real_8*REAL(nze_global, real_8)/ & (REAL(nfullrows, real_8)*REAL(nfullrows + 1, real_8) + & SUM(REAL(row_blk_size, real_8)*REAL(row_blk_size - 1, real_8))) END IF ELSE IF (nze_global .EQ. INT(nfullrows, KIND=int_8)*INT(nfullcols, KIND=int_8)) THEN occupation = 1.0_real_8 ELSE occupation = REAL(nze_global, real_8)/(REAL(nfullrows, real_8)*REAL(nfullcols, real_8)) END IF END IF ELSE occupation = 0.0_real_8 END IF END FUNCTION dbcsr_get_occupation SUBROUTINE dbcsr_clear(matrix) !! Clear a matrix (remove all blocks) TYPE(dbcsr_type), INTENT(INOUT) :: matrix TYPE(dbcsr_type) :: matrix_tmp CALL dbcsr_create(matrix_tmp, matrix) CALL dbcsr_release(matrix) matrix = matrix_tmp END SUBROUTINE SUBROUTINE dbcsr_trace_sd(matrix_a, trace) !! Trace of DBCSR matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! DBCSR matrices REAL(kind=real_8), INTENT(INOUT) :: trace !! the trace of the product of the matrices CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_sd' INTEGER :: handle REAL(kind=real_4) :: trace_4 CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8) THEN CALL dbcsr_trace_d(matrix_a, trace) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4) THEN trace_4 = 0.0_real_4 CALL dbcsr_trace_s(matrix_a, trace_4) trace = REAL(trace_4, real_8) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_trace_sd SUBROUTINE dbcsr_dot_sd(matrix_a, matrix_b, trace) !! Dot product of DBCSR matrices !! \result the dot product of the matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices REAL(kind=real_8), INTENT(INOUT) :: trace CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_dot_sd' INTEGER :: handle REAL(kind=real_4) :: trace_4 CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_8 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_8) THEN CALL dbcsr_dot_d(matrix_a, matrix_b, trace) ELSEIF (dbcsr_get_data_type(matrix_a) .EQ. dbcsr_type_real_4 .AND. & dbcsr_get_data_type(matrix_b) .EQ. dbcsr_type_real_4) THEN trace_4 = 0.0_real_4 CALL dbcsr_dot_s(matrix_a, matrix_b, trace_4) trace = REAL(trace_4, real_8) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") END IF CALL timestop(handle) END SUBROUTINE dbcsr_dot_sd # 1 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" 1 # 9 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 11 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 169 "/__w/dbcsr/dbcsr/src/ops/../data/dbcsr.fypp" # 2638 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" 2 # 2639 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" SUBROUTINE dbcsr_trace_d (matrix_a, trace) !! traces a DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! DBCSR matrix REAL(kind=real_8), INTENT(INOUT) :: trace !! the trace of the matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_d' INTEGER :: a_blk, a_col, a_col_size, & a_nze, a_row, a_row_size, i, & mynode, error_handle INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size, & row_dist, col_dist REAL(kind=real_8), DIMENSION(:), POINTER :: a_data, data_p INTEGER, DIMENSION(:, :), POINTER :: pgrid TYPE(dbcsr_distribution_obj) :: dist ! --------------------------------------------------------------------------- CALL timeset(routineN, error_handle) row_blk_size => array_data(matrix_a%row_blk_size) col_blk_size => array_data(matrix_a%col_blk_size) IF (dbcsr_get_data_type(matrix_a) /= dbcsr_type_real_8) & DBCSR_ABORT("Incompatible data types") CALL dbcsr_get_data(matrix_a%data_area, data_p) dist = dbcsr_distribution(matrix_a) mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dist)) pgrid => dbcsr_mp_pgrid(dbcsr_distribution_mp(dist)) row_dist => dbcsr_distribution_row_dist(dist) col_dist => dbcsr_distribution_col_dist(dist) ! ! let's go trace = REAL(0.0, real_8) DO a_row = 1, matrix_a%nblkrows_total a_row_size = row_blk_size(a_row) DO a_blk = matrix_a%row_p(a_row) + 1, matrix_a%row_p(a_row + 1) IF (a_blk .EQ. 0) CYCLE a_col = matrix_a%col_i(a_blk) IF (a_col .ne. a_row) CYCLE ! We must skip non-local blocks in a replicated matrix. IF (matrix_a%replication_type .NE. dbcsr_repl_full) THEN IF (mynode .NE. checker_square_proc(a_row, a_col, pgrid, row_dist, col_dist)) & CYCLE END IF a_col_size = col_blk_size(a_col) IF (a_row_size .NE. a_col_size) & DBCSR_ABORT("is that a square matrix?") a_nze = a_row_size**2 a_data => pointer_view(data_p, ABS(matrix_a%blk_p(a_blk)), & ABS(matrix_a%blk_p(a_blk)) + a_nze - 1) !data_a => matrix_a%data(ABS(matrix_a%blk_p(a_blk)):ABS(matrix_a%blk_p(a_blk))+a_nze-1) ! ! let's trace the block DO i = 1, a_row_size trace = trace + a_data((i - 1)*a_row_size + i) END DO END DO ! a_col END DO ! a_row ! ! summe CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) CALL timestop(error_handle) END SUBROUTINE dbcsr_trace_d SUBROUTINE dbcsr_dot_d (matrix_a, matrix_b, trace) !! Dot product of DBCSR matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices REAL(kind=real_8), INTENT(INOUT) :: trace !! the trace of the product of the matrices INTEGER :: a_blk, a_col, a_col_size, a_row_size, b_blk, b_col_size, & b_frst_blk, b_last_blk, b_row_size, nze, row, a_beg, a_end, b_beg, b_end CHARACTER :: matrix_a_type, matrix_b_type INTEGER, DIMENSION(:), POINTER :: a_col_blk_size, & a_row_blk_size, & b_col_blk_size, b_row_blk_size REAL(kind=real_8) :: sym_fac, fac LOGICAL :: found, matrix_a_symm, matrix_b_symm REAL(kind=real_8), DIMENSION(:), POINTER :: a_data, b_data ! --------------------------------------------------------------------------- IF (matrix_a%replication_type .NE. dbcsr_repl_none & .OR. matrix_b%replication_type .NE. dbcsr_repl_none) & DBCSR_ABORT("Trace of product of replicated matrices not yet possible.") sym_fac = REAL(1.0, real_8) matrix_a_type = dbcsr_get_matrix_type(matrix_a) matrix_b_type = dbcsr_get_matrix_type(matrix_b) matrix_a_symm = matrix_a_type == dbcsr_type_symmetric .OR. matrix_a_type == dbcsr_type_antisymmetric matrix_b_symm = matrix_b_type == dbcsr_type_symmetric .OR. matrix_b_type == dbcsr_type_antisymmetric IF (matrix_a_symm .AND. matrix_b_symm) sym_fac = REAL(2.0, real_8) ! tracing a symmetric with a general matrix is not implemented, as it would require communication of blocks IF (matrix_a_symm .NEQV. matrix_b_symm) & DBCSR_ABORT("Tracing general with symmetric matrix NYI") a_row_blk_size => array_data(matrix_a%row_blk_size) a_col_blk_size => array_data(matrix_a%col_blk_size) b_row_blk_size => array_data(matrix_b%row_blk_size) b_col_blk_size => array_data(matrix_b%col_blk_size) CALL dbcsr_get_data(matrix_a%data_area, a_data) CALL dbcsr_get_data(matrix_b%data_area, b_data) ! let's go trace = REAL(0.0, real_8) IF (matrix_a%nblkrows_total .NE. matrix_b%nblkrows_total) & DBCSR_ABORT("this combination of transpose is NYI") DO row = 1, matrix_a%nblkrows_total a_row_size = a_row_blk_size(row) b_row_size = b_row_blk_size(row) IF (a_row_size .NE. b_row_size) DBCSR_ABORT("matrices not consistent") b_blk = matrix_b%row_p(row) + 1 b_frst_blk = matrix_b%row_p(row) + 1 b_last_blk = matrix_b%row_p(row + 1) DO a_blk = matrix_a%row_p(row) + 1, matrix_a%row_p(row + 1) IF (matrix_a%blk_p(a_blk) .EQ. 0) CYCLE ! Deleted block a_col = matrix_a%col_i(a_blk) a_col_size = a_col_blk_size(a_col) ! ! find the b_blk we assume here that the columns are ordered ! CALL dbcsr_find_column(a_col, b_frst_blk, b_last_blk, matrix_b%col_i, & matrix_b%blk_p, b_blk, found) IF (found) THEN b_col_size = b_col_blk_size(a_col) IF (a_col_size .NE. b_col_size) DBCSR_ABORT("matrices not consistent") ! nze = a_row_size*a_col_size ! IF (nze .GT. 0) THEN ! ! let's trace the blocks a_beg = ABS(matrix_a%blk_p(a_blk)) a_end = a_beg + nze - 1 b_beg = ABS(matrix_b%blk_p(b_blk)) b_end = b_beg + nze - 1 fac = REAL(1.0, real_8) IF (row .NE. a_col) fac = sym_fac trace = trace + fac*SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) END IF END IF END DO ! a_col END DO ! a_row ! ! sum CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) END SUBROUTINE dbcsr_dot_d SUBROUTINE dbcsr_scale_d (matrix_a, alpha_scalar, last_column) !! Interface for matrix scaling by a scalar TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a REAL(kind=real_8), INTENT(IN) :: alpha_scalar INTEGER, INTENT(IN), OPTIONAL :: last_column CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_d' INTEGER :: error_handler TYPE(dbcsr_scalar_type) :: sc sc = dbcsr_scalar(alpha_scalar) CALL dbcsr_scalar_fill_all(sc) sc%data_type = dbcsr_get_data_type(matrix_a) CALL timeset(routineN, error_handler) IF (PRESENT(last_column)) THEN CALL dbcsr_scale_anytype(matrix_a, & alpha_scalar=sc, & limits=(/0, 0, 0, last_column/)) ELSE CALL dbcsr_scale_anytype(matrix_a, alpha_scalar=sc) END IF CALL timestop(error_handler) END SUBROUTINE dbcsr_scale_d SUBROUTINE dbcsr_scale_by_vector_d (matrix_a, alpha, side) !! Interface for matrix scaling by a vector TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a REAL(kind=real_8), DIMENSION(:), INTENT(IN), TARGET, CONTIGUOUS :: alpha CHARACTER(LEN=*), INTENT(IN) :: side REAL(kind=real_8), DIMENSION(:), POINTER, CONTIGUOUS :: tmp_p TYPE(dbcsr_data_obj) :: enc_alpha_vec CALL dbcsr_data_init(enc_alpha_vec) CALL dbcsr_data_new(enc_alpha_vec, dbcsr_type_real_8) tmp_p => alpha CALL dbcsr_data_set_pointer(enc_alpha_vec, tmp_p) CALL dbcsr_scale_by_vector_anytype(matrix_a, enc_alpha_vec, side) CALL dbcsr_data_clear_pointer(enc_alpha_vec) CALL dbcsr_data_release(enc_alpha_vec) END SUBROUTINE dbcsr_scale_by_vector_d SUBROUTINE dbcsr_set_d (matrix, alpha) !! Interface for dbcsr_set TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_8), INTENT(IN) :: alpha CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set' INTEGER :: col, handle, row TYPE(dbcsr_iterator) :: iter REAL(kind=real_8), DIMENSION(:, :), POINTER :: block LOGICAL :: tr CALL timeset(routineN, handle) IF (alpha == 0.0_real_8) THEN CALL dbcsr_zero(matrix) ELSE IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_8) & DBCSR_ABORT("Incompatible data types") !TODO: could be speedup by direct assignment to data_area, similar to dbcsr_zero() CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, block, tr) block(:, :) = alpha END DO CALL dbcsr_iterator_stop(iter) END IF CALL timestop(handle) END SUBROUTINE dbcsr_set_d SUBROUTINE dbcsr_filter_d (matrix, eps, method, use_absolute, & filter_diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_8), INTENT(IN) :: eps INTEGER, INTENT(IN), OPTIONAL :: method LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag CALL dbcsr_filter_anytype(matrix, dbcsr_scalar(eps), method, & use_absolute, filter_diag) END SUBROUTINE dbcsr_filter_d SUBROUTINE dbcsr_set_diag_d (matrix, diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_8), DIMENSION(:), INTENT(IN) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter REAL(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_8) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) block(i, i) = diag(row_offset + i - 1) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_set_diag_d SUBROUTINE dbcsr_get_diag_d (matrix, diag) TYPE(dbcsr_type), INTENT(IN) :: matrix REAL(kind=real_8), DIMENSION(:), INTENT(OUT) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter REAL(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_8) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") diag(:) = 0.0_real_8 CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) diag(row_offset + i - 1) = block(i, i) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_get_diag_d SUBROUTINE dbcsr_add_on_diag_d (matrix, alpha) !! add a constant to the diagonal of a matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix REAL(kind=real_8), INTENT(IN) :: alpha !! scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag' INTEGER :: handle, mynode, node, irow, i, row_size LOGICAL :: found, tr REAL(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_8) & DBCSR_ABORT("Incompatible data types") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) DO irow = 1, dbcsr_nblkrows_total(matrix) CALL dbcsr_get_stored_coordinates(matrix, irow, irow, node) IF (node /= mynode) CYCLE CALL dbcsr_get_block_p(matrix, irow, irow, block, tr, found, row_size=row_size) IF (.NOT. found) THEN ALLOCATE (block(row_size, row_size)) block(:, :) = 0.0_real_8 END IF DO i = 1, row_size block(i, i) = block(i, i) + alpha END DO IF (.NOT. found) THEN CALL dbcsr_put_block(matrix, irow, irow, block) DEALLOCATE (block) END IF END DO CALL dbcsr_finalize(matrix) CALL timestop(handle) END SUBROUTINE dbcsr_add_on_diag_d SUBROUTINE dbcsr_update_contiguous_blocks_d (matrix_a, matrix_b, first_lb_a, first_lb_b, nze, & do_scale, my_beta_scalar, found, iw) !! Low level function to sum contiguous chunks of blocks of the matrices (matrix_a = matrix_a + beta*matrix_b) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER, INTENT(IN) :: first_lb_a, first_lb_b, nze, iw LOGICAL, INTENT(IN) :: found, do_scale INTEGER :: ub_a, ub_b ub_a = first_lb_a + nze - 1 ub_b = first_lb_b + nze - 1 IF (found) THEN IF (do_scale) THEN CALL daxpy(nze, my_beta_scalar%r_dp, & matrix_b%data_area%d%r_dp (first_lb_b:ub_b), 1, & matrix_a%data_area%d%r_dp (first_lb_a:ub_a), 1) ELSE matrix_a%data_area%d%r_dp (first_lb_a:ub_a) = & matrix_a%data_area%d%r_dp (first_lb_a:ub_a) + & matrix_b%data_area%d%r_dp (first_lb_b:ub_b) END IF ELSE IF (do_scale) THEN matrix_a%wms(iw)%data_area%d%r_dp (first_lb_a:ub_a) = & my_beta_scalar%r_dp* & matrix_b%data_area%d%r_dp (first_lb_b:ub_b) ELSE matrix_a%wms(iw)%data_area%d%r_dp (first_lb_a:ub_a) = & matrix_b%data_area%d%r_dp (first_lb_b:ub_b) END IF END IF END SUBROUTINE dbcsr_update_contiguous_blocks_d SUBROUTINE dbcsr_add_anytype_d (matrix_a, matrix_b, iter, iw, do_scale, & my_beta_scalar, my_flop) !! Low level function to sum two matrices (matrix_a = matrix_a + beta*matrix_b TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_iterator), INTENT(INOUT) :: iter INTEGER, INTENT(IN) :: iw LOGICAL, INTENT(IN) :: do_scale TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER(KIND=int_8), INTENT(INOUT) :: my_flop INTEGER :: row, col, row_size, col_size, & nze, tot_nze, blk, & lb_a, first_lb_a, lb_a_val, & lb_b, first_lb_b INTEGER, DIMENSION(2) :: lb_row_blk LOGICAL :: was_found, found, tr ! some start values lb_row_blk(:) = 0 first_lb_a = matrix_a%wms(iw)%datasize + 1 first_lb_b = 0 tot_nze = 0 ! DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, blk, tr, lb_b, row_size, col_size) nze = row_size*col_size IF (nze .LE. 0) CYCLE IF (lb_row_blk(1) .LT. row) THEN lb_row_blk(1) = row lb_row_blk(2) = matrix_a%row_p(row) + 1 END IF ! get b-block index lb_b = ABS(lb_b) CALL dbcsr_find_column(col, lb_row_blk(2), matrix_a%row_p(row + 1), matrix_a%col_i, matrix_a%blk_p, blk, found) lb_row_blk(2) = blk + 1 ! get index of a-block lb_a whether found (from matrix_a) or not (from workspace array) IF (found) THEN my_flop = my_flop + nze*2 lb_a = ABS(matrix_a%blk_p(blk)) ELSE lb_a = matrix_a%wms(iw)%datasize + 1 lb_a_val = lb_a IF (tr) lb_a_val = -lb_a matrix_a%wms(iw)%lastblk = matrix_a%wms(iw)%lastblk + 1 matrix_a%wms(iw)%row_i(matrix_a%wms(iw)%lastblk) = row matrix_a%wms(iw)%col_i(matrix_a%wms(iw)%lastblk) = col matrix_a%wms(iw)%blk_p(matrix_a%wms(iw)%lastblk) = lb_a_val matrix_a%wms(iw)%datasize = matrix_a%wms(iw)%datasize + nze END IF ! at the first iteration we skip this and go directly to initialization after IF (first_lb_b .NE. 0) THEN ! if found status is the same as before then probably we are in contiguous blocks IF ((found .EQV. was_found) .AND. & (first_lb_b + tot_nze .EQ. lb_b) .AND. & (first_lb_a + tot_nze) .EQ. lb_a) THEN tot_nze = tot_nze + nze CYCLE END IF ! save block chunk CALL dbcsr_update_contiguous_blocks_d (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF ! first_lb_a = lb_a first_lb_b = lb_b tot_nze = nze was_found = found END DO ! save the last block or chunk of blocks IF (first_lb_b .NE. 0) THEN call dbcsr_update_contiguous_blocks_d (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF END SUBROUTINE dbcsr_add_anytype_d # 2639 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" SUBROUTINE dbcsr_trace_s (matrix_a, trace) !! traces a DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! DBCSR matrix REAL(kind=real_4), INTENT(INOUT) :: trace !! the trace of the matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_s' INTEGER :: a_blk, a_col, a_col_size, & a_nze, a_row, a_row_size, i, & mynode, error_handle INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size, & row_dist, col_dist REAL(kind=real_4), DIMENSION(:), POINTER :: a_data, data_p INTEGER, DIMENSION(:, :), POINTER :: pgrid TYPE(dbcsr_distribution_obj) :: dist ! --------------------------------------------------------------------------- CALL timeset(routineN, error_handle) row_blk_size => array_data(matrix_a%row_blk_size) col_blk_size => array_data(matrix_a%col_blk_size) IF (dbcsr_get_data_type(matrix_a) /= dbcsr_type_real_4) & DBCSR_ABORT("Incompatible data types") CALL dbcsr_get_data(matrix_a%data_area, data_p) dist = dbcsr_distribution(matrix_a) mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dist)) pgrid => dbcsr_mp_pgrid(dbcsr_distribution_mp(dist)) row_dist => dbcsr_distribution_row_dist(dist) col_dist => dbcsr_distribution_col_dist(dist) ! ! let's go trace = REAL(0.0, real_4) DO a_row = 1, matrix_a%nblkrows_total a_row_size = row_blk_size(a_row) DO a_blk = matrix_a%row_p(a_row) + 1, matrix_a%row_p(a_row + 1) IF (a_blk .EQ. 0) CYCLE a_col = matrix_a%col_i(a_blk) IF (a_col .ne. a_row) CYCLE ! We must skip non-local blocks in a replicated matrix. IF (matrix_a%replication_type .NE. dbcsr_repl_full) THEN IF (mynode .NE. checker_square_proc(a_row, a_col, pgrid, row_dist, col_dist)) & CYCLE END IF a_col_size = col_blk_size(a_col) IF (a_row_size .NE. a_col_size) & DBCSR_ABORT("is that a square matrix?") a_nze = a_row_size**2 a_data => pointer_view(data_p, ABS(matrix_a%blk_p(a_blk)), & ABS(matrix_a%blk_p(a_blk)) + a_nze - 1) !data_a => matrix_a%data(ABS(matrix_a%blk_p(a_blk)):ABS(matrix_a%blk_p(a_blk))+a_nze-1) ! ! let's trace the block DO i = 1, a_row_size trace = trace + a_data((i - 1)*a_row_size + i) END DO END DO ! a_col END DO ! a_row ! ! summe CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) CALL timestop(error_handle) END SUBROUTINE dbcsr_trace_s SUBROUTINE dbcsr_dot_s (matrix_a, matrix_b, trace) !! Dot product of DBCSR matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices REAL(kind=real_4), INTENT(INOUT) :: trace !! the trace of the product of the matrices INTEGER :: a_blk, a_col, a_col_size, a_row_size, b_blk, b_col_size, & b_frst_blk, b_last_blk, b_row_size, nze, row, a_beg, a_end, b_beg, b_end CHARACTER :: matrix_a_type, matrix_b_type INTEGER, DIMENSION(:), POINTER :: a_col_blk_size, & a_row_blk_size, & b_col_blk_size, b_row_blk_size REAL(kind=real_4) :: sym_fac, fac LOGICAL :: found, matrix_a_symm, matrix_b_symm REAL(kind=real_4), DIMENSION(:), POINTER :: a_data, b_data ! --------------------------------------------------------------------------- IF (matrix_a%replication_type .NE. dbcsr_repl_none & .OR. matrix_b%replication_type .NE. dbcsr_repl_none) & DBCSR_ABORT("Trace of product of replicated matrices not yet possible.") sym_fac = REAL(1.0, real_4) matrix_a_type = dbcsr_get_matrix_type(matrix_a) matrix_b_type = dbcsr_get_matrix_type(matrix_b) matrix_a_symm = matrix_a_type == dbcsr_type_symmetric .OR. matrix_a_type == dbcsr_type_antisymmetric matrix_b_symm = matrix_b_type == dbcsr_type_symmetric .OR. matrix_b_type == dbcsr_type_antisymmetric IF (matrix_a_symm .AND. matrix_b_symm) sym_fac = REAL(2.0, real_4) ! tracing a symmetric with a general matrix is not implemented, as it would require communication of blocks IF (matrix_a_symm .NEQV. matrix_b_symm) & DBCSR_ABORT("Tracing general with symmetric matrix NYI") a_row_blk_size => array_data(matrix_a%row_blk_size) a_col_blk_size => array_data(matrix_a%col_blk_size) b_row_blk_size => array_data(matrix_b%row_blk_size) b_col_blk_size => array_data(matrix_b%col_blk_size) CALL dbcsr_get_data(matrix_a%data_area, a_data) CALL dbcsr_get_data(matrix_b%data_area, b_data) ! let's go trace = REAL(0.0, real_4) IF (matrix_a%nblkrows_total .NE. matrix_b%nblkrows_total) & DBCSR_ABORT("this combination of transpose is NYI") DO row = 1, matrix_a%nblkrows_total a_row_size = a_row_blk_size(row) b_row_size = b_row_blk_size(row) IF (a_row_size .NE. b_row_size) DBCSR_ABORT("matrices not consistent") b_blk = matrix_b%row_p(row) + 1 b_frst_blk = matrix_b%row_p(row) + 1 b_last_blk = matrix_b%row_p(row + 1) DO a_blk = matrix_a%row_p(row) + 1, matrix_a%row_p(row + 1) IF (matrix_a%blk_p(a_blk) .EQ. 0) CYCLE ! Deleted block a_col = matrix_a%col_i(a_blk) a_col_size = a_col_blk_size(a_col) ! ! find the b_blk we assume here that the columns are ordered ! CALL dbcsr_find_column(a_col, b_frst_blk, b_last_blk, matrix_b%col_i, & matrix_b%blk_p, b_blk, found) IF (found) THEN b_col_size = b_col_blk_size(a_col) IF (a_col_size .NE. b_col_size) DBCSR_ABORT("matrices not consistent") ! nze = a_row_size*a_col_size ! IF (nze .GT. 0) THEN ! ! let's trace the blocks a_beg = ABS(matrix_a%blk_p(a_blk)) a_end = a_beg + nze - 1 b_beg = ABS(matrix_b%blk_p(b_blk)) b_end = b_beg + nze - 1 fac = REAL(1.0, real_4) IF (row .NE. a_col) fac = sym_fac trace = trace + fac*SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) END IF END IF END DO ! a_col END DO ! a_row ! ! sum CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) END SUBROUTINE dbcsr_dot_s SUBROUTINE dbcsr_scale_s (matrix_a, alpha_scalar, last_column) !! Interface for matrix scaling by a scalar TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a REAL(kind=real_4), INTENT(IN) :: alpha_scalar INTEGER, INTENT(IN), OPTIONAL :: last_column CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_s' INTEGER :: error_handler TYPE(dbcsr_scalar_type) :: sc sc = dbcsr_scalar(alpha_scalar) CALL dbcsr_scalar_fill_all(sc) sc%data_type = dbcsr_get_data_type(matrix_a) CALL timeset(routineN, error_handler) IF (PRESENT(last_column)) THEN CALL dbcsr_scale_anytype(matrix_a, & alpha_scalar=sc, & limits=(/0, 0, 0, last_column/)) ELSE CALL dbcsr_scale_anytype(matrix_a, alpha_scalar=sc) END IF CALL timestop(error_handler) END SUBROUTINE dbcsr_scale_s SUBROUTINE dbcsr_scale_by_vector_s (matrix_a, alpha, side) !! Interface for matrix scaling by a vector TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a REAL(kind=real_4), DIMENSION(:), INTENT(IN), TARGET, CONTIGUOUS :: alpha CHARACTER(LEN=*), INTENT(IN) :: side REAL(kind=real_4), DIMENSION(:), POINTER, CONTIGUOUS :: tmp_p TYPE(dbcsr_data_obj) :: enc_alpha_vec CALL dbcsr_data_init(enc_alpha_vec) CALL dbcsr_data_new(enc_alpha_vec, dbcsr_type_real_4) tmp_p => alpha CALL dbcsr_data_set_pointer(enc_alpha_vec, tmp_p) CALL dbcsr_scale_by_vector_anytype(matrix_a, enc_alpha_vec, side) CALL dbcsr_data_clear_pointer(enc_alpha_vec) CALL dbcsr_data_release(enc_alpha_vec) END SUBROUTINE dbcsr_scale_by_vector_s SUBROUTINE dbcsr_set_s (matrix, alpha) !! Interface for dbcsr_set TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_4), INTENT(IN) :: alpha CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set' INTEGER :: col, handle, row TYPE(dbcsr_iterator) :: iter REAL(kind=real_4), DIMENSION(:, :), POINTER :: block LOGICAL :: tr CALL timeset(routineN, handle) IF (alpha == 0.0_real_4) THEN CALL dbcsr_zero(matrix) ELSE IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_4) & DBCSR_ABORT("Incompatible data types") !TODO: could be speedup by direct assignment to data_area, similar to dbcsr_zero() CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, block, tr) block(:, :) = alpha END DO CALL dbcsr_iterator_stop(iter) END IF CALL timestop(handle) END SUBROUTINE dbcsr_set_s SUBROUTINE dbcsr_filter_s (matrix, eps, method, use_absolute, & filter_diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_4), INTENT(IN) :: eps INTEGER, INTENT(IN), OPTIONAL :: method LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag CALL dbcsr_filter_anytype(matrix, dbcsr_scalar(eps), method, & use_absolute, filter_diag) END SUBROUTINE dbcsr_filter_s SUBROUTINE dbcsr_set_diag_s (matrix, diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix REAL(kind=real_4), DIMENSION(:), INTENT(IN) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter REAL(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_4) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) block(i, i) = diag(row_offset + i - 1) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_set_diag_s SUBROUTINE dbcsr_get_diag_s (matrix, diag) TYPE(dbcsr_type), INTENT(IN) :: matrix REAL(kind=real_4), DIMENSION(:), INTENT(OUT) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter REAL(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_4) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") diag(:) = 0.0_real_4 CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) diag(row_offset + i - 1) = block(i, i) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_get_diag_s SUBROUTINE dbcsr_add_on_diag_s (matrix, alpha) !! add a constant to the diagonal of a matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix REAL(kind=real_4), INTENT(IN) :: alpha !! scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag' INTEGER :: handle, mynode, node, irow, i, row_size LOGICAL :: found, tr REAL(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_real_4) & DBCSR_ABORT("Incompatible data types") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) DO irow = 1, dbcsr_nblkrows_total(matrix) CALL dbcsr_get_stored_coordinates(matrix, irow, irow, node) IF (node /= mynode) CYCLE CALL dbcsr_get_block_p(matrix, irow, irow, block, tr, found, row_size=row_size) IF (.NOT. found) THEN ALLOCATE (block(row_size, row_size)) block(:, :) = 0.0_real_4 END IF DO i = 1, row_size block(i, i) = block(i, i) + alpha END DO IF (.NOT. found) THEN CALL dbcsr_put_block(matrix, irow, irow, block) DEALLOCATE (block) END IF END DO CALL dbcsr_finalize(matrix) CALL timestop(handle) END SUBROUTINE dbcsr_add_on_diag_s SUBROUTINE dbcsr_update_contiguous_blocks_s (matrix_a, matrix_b, first_lb_a, first_lb_b, nze, & do_scale, my_beta_scalar, found, iw) !! Low level function to sum contiguous chunks of blocks of the matrices (matrix_a = matrix_a + beta*matrix_b) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER, INTENT(IN) :: first_lb_a, first_lb_b, nze, iw LOGICAL, INTENT(IN) :: found, do_scale INTEGER :: ub_a, ub_b ub_a = first_lb_a + nze - 1 ub_b = first_lb_b + nze - 1 IF (found) THEN IF (do_scale) THEN CALL saxpy(nze, my_beta_scalar%r_sp, & matrix_b%data_area%d%r_sp (first_lb_b:ub_b), 1, & matrix_a%data_area%d%r_sp (first_lb_a:ub_a), 1) ELSE matrix_a%data_area%d%r_sp (first_lb_a:ub_a) = & matrix_a%data_area%d%r_sp (first_lb_a:ub_a) + & matrix_b%data_area%d%r_sp (first_lb_b:ub_b) END IF ELSE IF (do_scale) THEN matrix_a%wms(iw)%data_area%d%r_sp (first_lb_a:ub_a) = & my_beta_scalar%r_sp* & matrix_b%data_area%d%r_sp (first_lb_b:ub_b) ELSE matrix_a%wms(iw)%data_area%d%r_sp (first_lb_a:ub_a) = & matrix_b%data_area%d%r_sp (first_lb_b:ub_b) END IF END IF END SUBROUTINE dbcsr_update_contiguous_blocks_s SUBROUTINE dbcsr_add_anytype_s (matrix_a, matrix_b, iter, iw, do_scale, & my_beta_scalar, my_flop) !! Low level function to sum two matrices (matrix_a = matrix_a + beta*matrix_b TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_iterator), INTENT(INOUT) :: iter INTEGER, INTENT(IN) :: iw LOGICAL, INTENT(IN) :: do_scale TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER(KIND=int_8), INTENT(INOUT) :: my_flop INTEGER :: row, col, row_size, col_size, & nze, tot_nze, blk, & lb_a, first_lb_a, lb_a_val, & lb_b, first_lb_b INTEGER, DIMENSION(2) :: lb_row_blk LOGICAL :: was_found, found, tr ! some start values lb_row_blk(:) = 0 first_lb_a = matrix_a%wms(iw)%datasize + 1 first_lb_b = 0 tot_nze = 0 ! DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, blk, tr, lb_b, row_size, col_size) nze = row_size*col_size IF (nze .LE. 0) CYCLE IF (lb_row_blk(1) .LT. row) THEN lb_row_blk(1) = row lb_row_blk(2) = matrix_a%row_p(row) + 1 END IF ! get b-block index lb_b = ABS(lb_b) CALL dbcsr_find_column(col, lb_row_blk(2), matrix_a%row_p(row + 1), matrix_a%col_i, matrix_a%blk_p, blk, found) lb_row_blk(2) = blk + 1 ! get index of a-block lb_a whether found (from matrix_a) or not (from workspace array) IF (found) THEN my_flop = my_flop + nze*2 lb_a = ABS(matrix_a%blk_p(blk)) ELSE lb_a = matrix_a%wms(iw)%datasize + 1 lb_a_val = lb_a IF (tr) lb_a_val = -lb_a matrix_a%wms(iw)%lastblk = matrix_a%wms(iw)%lastblk + 1 matrix_a%wms(iw)%row_i(matrix_a%wms(iw)%lastblk) = row matrix_a%wms(iw)%col_i(matrix_a%wms(iw)%lastblk) = col matrix_a%wms(iw)%blk_p(matrix_a%wms(iw)%lastblk) = lb_a_val matrix_a%wms(iw)%datasize = matrix_a%wms(iw)%datasize + nze END IF ! at the first iteration we skip this and go directly to initialization after IF (first_lb_b .NE. 0) THEN ! if found status is the same as before then probably we are in contiguous blocks IF ((found .EQV. was_found) .AND. & (first_lb_b + tot_nze .EQ. lb_b) .AND. & (first_lb_a + tot_nze) .EQ. lb_a) THEN tot_nze = tot_nze + nze CYCLE END IF ! save block chunk CALL dbcsr_update_contiguous_blocks_s (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF ! first_lb_a = lb_a first_lb_b = lb_b tot_nze = nze was_found = found END DO ! save the last block or chunk of blocks IF (first_lb_b .NE. 0) THEN call dbcsr_update_contiguous_blocks_s (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF END SUBROUTINE dbcsr_add_anytype_s # 2639 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" SUBROUTINE dbcsr_trace_z (matrix_a, trace) !! traces a DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! DBCSR matrix COMPLEX(kind=real_8), INTENT(INOUT) :: trace !! the trace of the matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_z' INTEGER :: a_blk, a_col, a_col_size, & a_nze, a_row, a_row_size, i, & mynode, error_handle INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size, & row_dist, col_dist COMPLEX(kind=real_8), DIMENSION(:), POINTER :: a_data, data_p INTEGER, DIMENSION(:, :), POINTER :: pgrid TYPE(dbcsr_distribution_obj) :: dist ! --------------------------------------------------------------------------- CALL timeset(routineN, error_handle) row_blk_size => array_data(matrix_a%row_blk_size) col_blk_size => array_data(matrix_a%col_blk_size) IF (dbcsr_get_data_type(matrix_a) /= dbcsr_type_complex_8) & DBCSR_ABORT("Incompatible data types") CALL dbcsr_get_data(matrix_a%data_area, data_p) dist = dbcsr_distribution(matrix_a) mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dist)) pgrid => dbcsr_mp_pgrid(dbcsr_distribution_mp(dist)) row_dist => dbcsr_distribution_row_dist(dist) col_dist => dbcsr_distribution_col_dist(dist) ! ! let's go trace = REAL(0.0, real_8) DO a_row = 1, matrix_a%nblkrows_total a_row_size = row_blk_size(a_row) DO a_blk = matrix_a%row_p(a_row) + 1, matrix_a%row_p(a_row + 1) IF (a_blk .EQ. 0) CYCLE a_col = matrix_a%col_i(a_blk) IF (a_col .ne. a_row) CYCLE ! We must skip non-local blocks in a replicated matrix. IF (matrix_a%replication_type .NE. dbcsr_repl_full) THEN IF (mynode .NE. checker_square_proc(a_row, a_col, pgrid, row_dist, col_dist)) & CYCLE END IF a_col_size = col_blk_size(a_col) IF (a_row_size .NE. a_col_size) & DBCSR_ABORT("is that a square matrix?") a_nze = a_row_size**2 a_data => pointer_view(data_p, ABS(matrix_a%blk_p(a_blk)), & ABS(matrix_a%blk_p(a_blk)) + a_nze - 1) !data_a => matrix_a%data(ABS(matrix_a%blk_p(a_blk)):ABS(matrix_a%blk_p(a_blk))+a_nze-1) ! ! let's trace the block DO i = 1, a_row_size trace = trace + a_data((i - 1)*a_row_size + i) END DO END DO ! a_col END DO ! a_row ! ! summe CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) CALL timestop(error_handle) END SUBROUTINE dbcsr_trace_z SUBROUTINE dbcsr_dot_z (matrix_a, matrix_b, trace) !! Dot product of DBCSR matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices COMPLEX(kind=real_8), INTENT(INOUT) :: trace !! the trace of the product of the matrices INTEGER :: a_blk, a_col, a_col_size, a_row_size, b_blk, b_col_size, & b_frst_blk, b_last_blk, b_row_size, nze, row, a_beg, a_end, b_beg, b_end CHARACTER :: matrix_a_type, matrix_b_type INTEGER, DIMENSION(:), POINTER :: a_col_blk_size, & a_row_blk_size, & b_col_blk_size, b_row_blk_size COMPLEX(kind=real_8) :: sym_fac, fac LOGICAL :: found, matrix_a_symm, matrix_b_symm COMPLEX(kind=real_8), DIMENSION(:), POINTER :: a_data, b_data ! --------------------------------------------------------------------------- IF (matrix_a%replication_type .NE. dbcsr_repl_none & .OR. matrix_b%replication_type .NE. dbcsr_repl_none) & DBCSR_ABORT("Trace of product of replicated matrices not yet possible.") sym_fac = REAL(1.0, real_8) matrix_a_type = dbcsr_get_matrix_type(matrix_a) matrix_b_type = dbcsr_get_matrix_type(matrix_b) matrix_a_symm = matrix_a_type == dbcsr_type_symmetric .OR. matrix_a_type == dbcsr_type_antisymmetric matrix_b_symm = matrix_b_type == dbcsr_type_symmetric .OR. matrix_b_type == dbcsr_type_antisymmetric IF (matrix_a_symm .AND. matrix_b_symm) sym_fac = REAL(2.0, real_8) ! tracing a symmetric with a general matrix is not implemented, as it would require communication of blocks IF (matrix_a_symm .NEQV. matrix_b_symm) & DBCSR_ABORT("Tracing general with symmetric matrix NYI") a_row_blk_size => array_data(matrix_a%row_blk_size) a_col_blk_size => array_data(matrix_a%col_blk_size) b_row_blk_size => array_data(matrix_b%row_blk_size) b_col_blk_size => array_data(matrix_b%col_blk_size) CALL dbcsr_get_data(matrix_a%data_area, a_data) CALL dbcsr_get_data(matrix_b%data_area, b_data) ! let's go trace = REAL(0.0, real_8) IF (matrix_a%nblkrows_total .NE. matrix_b%nblkrows_total) & DBCSR_ABORT("this combination of transpose is NYI") DO row = 1, matrix_a%nblkrows_total a_row_size = a_row_blk_size(row) b_row_size = b_row_blk_size(row) IF (a_row_size .NE. b_row_size) DBCSR_ABORT("matrices not consistent") b_blk = matrix_b%row_p(row) + 1 b_frst_blk = matrix_b%row_p(row) + 1 b_last_blk = matrix_b%row_p(row + 1) DO a_blk = matrix_a%row_p(row) + 1, matrix_a%row_p(row + 1) IF (matrix_a%blk_p(a_blk) .EQ. 0) CYCLE ! Deleted block a_col = matrix_a%col_i(a_blk) a_col_size = a_col_blk_size(a_col) ! ! find the b_blk we assume here that the columns are ordered ! CALL dbcsr_find_column(a_col, b_frst_blk, b_last_blk, matrix_b%col_i, & matrix_b%blk_p, b_blk, found) IF (found) THEN b_col_size = b_col_blk_size(a_col) IF (a_col_size .NE. b_col_size) DBCSR_ABORT("matrices not consistent") ! nze = a_row_size*a_col_size ! IF (nze .GT. 0) THEN ! ! let's trace the blocks a_beg = ABS(matrix_a%blk_p(a_blk)) a_end = a_beg + nze - 1 b_beg = ABS(matrix_b%blk_p(b_blk)) b_end = b_beg + nze - 1 fac = REAL(1.0, real_8) IF (row .NE. a_col) fac = sym_fac trace = trace + fac*SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) END IF END IF END DO ! a_col END DO ! a_row ! ! sum CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) END SUBROUTINE dbcsr_dot_z SUBROUTINE dbcsr_scale_z (matrix_a, alpha_scalar, last_column) !! Interface for matrix scaling by a scalar TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a COMPLEX(kind=real_8), INTENT(IN) :: alpha_scalar INTEGER, INTENT(IN), OPTIONAL :: last_column CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_z' INTEGER :: error_handler TYPE(dbcsr_scalar_type) :: sc sc = dbcsr_scalar(alpha_scalar) CALL dbcsr_scalar_fill_all(sc) sc%data_type = dbcsr_get_data_type(matrix_a) CALL timeset(routineN, error_handler) IF (PRESENT(last_column)) THEN CALL dbcsr_scale_anytype(matrix_a, & alpha_scalar=sc, & limits=(/0, 0, 0, last_column/)) ELSE CALL dbcsr_scale_anytype(matrix_a, alpha_scalar=sc) END IF CALL timestop(error_handler) END SUBROUTINE dbcsr_scale_z SUBROUTINE dbcsr_scale_by_vector_z (matrix_a, alpha, side) !! Interface for matrix scaling by a vector TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a COMPLEX(kind=real_8), DIMENSION(:), INTENT(IN), TARGET, CONTIGUOUS :: alpha CHARACTER(LEN=*), INTENT(IN) :: side COMPLEX(kind=real_8), DIMENSION(:), POINTER, CONTIGUOUS :: tmp_p TYPE(dbcsr_data_obj) :: enc_alpha_vec CALL dbcsr_data_init(enc_alpha_vec) CALL dbcsr_data_new(enc_alpha_vec, dbcsr_type_complex_8) tmp_p => alpha CALL dbcsr_data_set_pointer(enc_alpha_vec, tmp_p) CALL dbcsr_scale_by_vector_anytype(matrix_a, enc_alpha_vec, side) CALL dbcsr_data_clear_pointer(enc_alpha_vec) CALL dbcsr_data_release(enc_alpha_vec) END SUBROUTINE dbcsr_scale_by_vector_z SUBROUTINE dbcsr_set_z (matrix, alpha) !! Interface for dbcsr_set TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_8), INTENT(IN) :: alpha CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set' INTEGER :: col, handle, row TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_8), DIMENSION(:, :), POINTER :: block LOGICAL :: tr CALL timeset(routineN, handle) IF (alpha == CMPLX(0.0, 0.0, real_8)) THEN CALL dbcsr_zero(matrix) ELSE IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_8) & DBCSR_ABORT("Incompatible data types") !TODO: could be speedup by direct assignment to data_area, similar to dbcsr_zero() CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, block, tr) block(:, :) = alpha END DO CALL dbcsr_iterator_stop(iter) END IF CALL timestop(handle) END SUBROUTINE dbcsr_set_z SUBROUTINE dbcsr_filter_z (matrix, eps, method, use_absolute, & filter_diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_8), INTENT(IN) :: eps INTEGER, INTENT(IN), OPTIONAL :: method LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag CALL dbcsr_filter_anytype(matrix, dbcsr_scalar(eps), method, & use_absolute, filter_diag) END SUBROUTINE dbcsr_filter_z SUBROUTINE dbcsr_set_diag_z (matrix, diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_8), DIMENSION(:), INTENT(IN) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_8) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) block(i, i) = diag(row_offset + i - 1) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_set_diag_z SUBROUTINE dbcsr_get_diag_z (matrix, diag) TYPE(dbcsr_type), INTENT(IN) :: matrix COMPLEX(kind=real_8), DIMENSION(:), INTENT(OUT) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_8) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") diag(:) = CMPLX(0.0, 0.0, real_8) CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) diag(row_offset + i - 1) = block(i, i) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_get_diag_z SUBROUTINE dbcsr_add_on_diag_z (matrix, alpha) !! add a constant to the diagonal of a matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix COMPLEX(kind=real_8), INTENT(IN) :: alpha !! scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag' INTEGER :: handle, mynode, node, irow, i, row_size LOGICAL :: found, tr COMPLEX(kind=real_8), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_8) & DBCSR_ABORT("Incompatible data types") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) DO irow = 1, dbcsr_nblkrows_total(matrix) CALL dbcsr_get_stored_coordinates(matrix, irow, irow, node) IF (node /= mynode) CYCLE CALL dbcsr_get_block_p(matrix, irow, irow, block, tr, found, row_size=row_size) IF (.NOT. found) THEN ALLOCATE (block(row_size, row_size)) block(:, :) = CMPLX(0.0, 0.0, real_8) END IF DO i = 1, row_size block(i, i) = block(i, i) + alpha END DO IF (.NOT. found) THEN CALL dbcsr_put_block(matrix, irow, irow, block) DEALLOCATE (block) END IF END DO CALL dbcsr_finalize(matrix) CALL timestop(handle) END SUBROUTINE dbcsr_add_on_diag_z SUBROUTINE dbcsr_update_contiguous_blocks_z (matrix_a, matrix_b, first_lb_a, first_lb_b, nze, & do_scale, my_beta_scalar, found, iw) !! Low level function to sum contiguous chunks of blocks of the matrices (matrix_a = matrix_a + beta*matrix_b) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER, INTENT(IN) :: first_lb_a, first_lb_b, nze, iw LOGICAL, INTENT(IN) :: found, do_scale INTEGER :: ub_a, ub_b ub_a = first_lb_a + nze - 1 ub_b = first_lb_b + nze - 1 IF (found) THEN IF (do_scale) THEN CALL zaxpy(nze, my_beta_scalar%c_dp, & matrix_b%data_area%d%c_dp (first_lb_b:ub_b), 1, & matrix_a%data_area%d%c_dp (first_lb_a:ub_a), 1) ELSE matrix_a%data_area%d%c_dp (first_lb_a:ub_a) = & matrix_a%data_area%d%c_dp (first_lb_a:ub_a) + & matrix_b%data_area%d%c_dp (first_lb_b:ub_b) END IF ELSE IF (do_scale) THEN matrix_a%wms(iw)%data_area%d%c_dp (first_lb_a:ub_a) = & my_beta_scalar%c_dp* & matrix_b%data_area%d%c_dp (first_lb_b:ub_b) ELSE matrix_a%wms(iw)%data_area%d%c_dp (first_lb_a:ub_a) = & matrix_b%data_area%d%c_dp (first_lb_b:ub_b) END IF END IF END SUBROUTINE dbcsr_update_contiguous_blocks_z SUBROUTINE dbcsr_add_anytype_z (matrix_a, matrix_b, iter, iw, do_scale, & my_beta_scalar, my_flop) !! Low level function to sum two matrices (matrix_a = matrix_a + beta*matrix_b TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_iterator), INTENT(INOUT) :: iter INTEGER, INTENT(IN) :: iw LOGICAL, INTENT(IN) :: do_scale TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER(KIND=int_8), INTENT(INOUT) :: my_flop INTEGER :: row, col, row_size, col_size, & nze, tot_nze, blk, & lb_a, first_lb_a, lb_a_val, & lb_b, first_lb_b INTEGER, DIMENSION(2) :: lb_row_blk LOGICAL :: was_found, found, tr ! some start values lb_row_blk(:) = 0 first_lb_a = matrix_a%wms(iw)%datasize + 1 first_lb_b = 0 tot_nze = 0 ! DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, blk, tr, lb_b, row_size, col_size) nze = row_size*col_size IF (nze .LE. 0) CYCLE IF (lb_row_blk(1) .LT. row) THEN lb_row_blk(1) = row lb_row_blk(2) = matrix_a%row_p(row) + 1 END IF ! get b-block index lb_b = ABS(lb_b) CALL dbcsr_find_column(col, lb_row_blk(2), matrix_a%row_p(row + 1), matrix_a%col_i, matrix_a%blk_p, blk, found) lb_row_blk(2) = blk + 1 ! get index of a-block lb_a whether found (from matrix_a) or not (from workspace array) IF (found) THEN my_flop = my_flop + nze*2 lb_a = ABS(matrix_a%blk_p(blk)) ELSE lb_a = matrix_a%wms(iw)%datasize + 1 lb_a_val = lb_a IF (tr) lb_a_val = -lb_a matrix_a%wms(iw)%lastblk = matrix_a%wms(iw)%lastblk + 1 matrix_a%wms(iw)%row_i(matrix_a%wms(iw)%lastblk) = row matrix_a%wms(iw)%col_i(matrix_a%wms(iw)%lastblk) = col matrix_a%wms(iw)%blk_p(matrix_a%wms(iw)%lastblk) = lb_a_val matrix_a%wms(iw)%datasize = matrix_a%wms(iw)%datasize + nze END IF ! at the first iteration we skip this and go directly to initialization after IF (first_lb_b .NE. 0) THEN ! if found status is the same as before then probably we are in contiguous blocks IF ((found .EQV. was_found) .AND. & (first_lb_b + tot_nze .EQ. lb_b) .AND. & (first_lb_a + tot_nze) .EQ. lb_a) THEN tot_nze = tot_nze + nze CYCLE END IF ! save block chunk CALL dbcsr_update_contiguous_blocks_z (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF ! first_lb_a = lb_a first_lb_b = lb_b tot_nze = nze was_found = found END DO ! save the last block or chunk of blocks IF (first_lb_b .NE. 0) THEN call dbcsr_update_contiguous_blocks_z (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF END SUBROUTINE dbcsr_add_anytype_z # 2639 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" SUBROUTINE dbcsr_trace_c (matrix_a, trace) !! traces a DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_a !! DBCSR matrix COMPLEX(kind=real_4), INTENT(INOUT) :: trace !! the trace of the matrix CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_trace_c' INTEGER :: a_blk, a_col, a_col_size, & a_nze, a_row, a_row_size, i, & mynode, error_handle INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size, & row_dist, col_dist COMPLEX(kind=real_4), DIMENSION(:), POINTER :: a_data, data_p INTEGER, DIMENSION(:, :), POINTER :: pgrid TYPE(dbcsr_distribution_obj) :: dist ! --------------------------------------------------------------------------- CALL timeset(routineN, error_handle) row_blk_size => array_data(matrix_a%row_blk_size) col_blk_size => array_data(matrix_a%col_blk_size) IF (dbcsr_get_data_type(matrix_a) /= dbcsr_type_complex_4) & DBCSR_ABORT("Incompatible data types") CALL dbcsr_get_data(matrix_a%data_area, data_p) dist = dbcsr_distribution(matrix_a) mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dist)) pgrid => dbcsr_mp_pgrid(dbcsr_distribution_mp(dist)) row_dist => dbcsr_distribution_row_dist(dist) col_dist => dbcsr_distribution_col_dist(dist) ! ! let's go trace = REAL(0.0, real_4) DO a_row = 1, matrix_a%nblkrows_total a_row_size = row_blk_size(a_row) DO a_blk = matrix_a%row_p(a_row) + 1, matrix_a%row_p(a_row + 1) IF (a_blk .EQ. 0) CYCLE a_col = matrix_a%col_i(a_blk) IF (a_col .ne. a_row) CYCLE ! We must skip non-local blocks in a replicated matrix. IF (matrix_a%replication_type .NE. dbcsr_repl_full) THEN IF (mynode .NE. checker_square_proc(a_row, a_col, pgrid, row_dist, col_dist)) & CYCLE END IF a_col_size = col_blk_size(a_col) IF (a_row_size .NE. a_col_size) & DBCSR_ABORT("is that a square matrix?") a_nze = a_row_size**2 a_data => pointer_view(data_p, ABS(matrix_a%blk_p(a_blk)), & ABS(matrix_a%blk_p(a_blk)) + a_nze - 1) !data_a => matrix_a%data(ABS(matrix_a%blk_p(a_blk)):ABS(matrix_a%blk_p(a_blk))+a_nze-1) ! ! let's trace the block DO i = 1, a_row_size trace = trace + a_data((i - 1)*a_row_size + i) END DO END DO ! a_col END DO ! a_row ! ! summe CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) CALL timestop(error_handle) END SUBROUTINE dbcsr_trace_c SUBROUTINE dbcsr_dot_c (matrix_a, matrix_b, trace) !! Dot product of DBCSR matrices TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices COMPLEX(kind=real_4), INTENT(INOUT) :: trace !! the trace of the product of the matrices INTEGER :: a_blk, a_col, a_col_size, a_row_size, b_blk, b_col_size, & b_frst_blk, b_last_blk, b_row_size, nze, row, a_beg, a_end, b_beg, b_end CHARACTER :: matrix_a_type, matrix_b_type INTEGER, DIMENSION(:), POINTER :: a_col_blk_size, & a_row_blk_size, & b_col_blk_size, b_row_blk_size COMPLEX(kind=real_4) :: sym_fac, fac LOGICAL :: found, matrix_a_symm, matrix_b_symm COMPLEX(kind=real_4), DIMENSION(:), POINTER :: a_data, b_data ! --------------------------------------------------------------------------- IF (matrix_a%replication_type .NE. dbcsr_repl_none & .OR. matrix_b%replication_type .NE. dbcsr_repl_none) & DBCSR_ABORT("Trace of product of replicated matrices not yet possible.") sym_fac = REAL(1.0, real_4) matrix_a_type = dbcsr_get_matrix_type(matrix_a) matrix_b_type = dbcsr_get_matrix_type(matrix_b) matrix_a_symm = matrix_a_type == dbcsr_type_symmetric .OR. matrix_a_type == dbcsr_type_antisymmetric matrix_b_symm = matrix_b_type == dbcsr_type_symmetric .OR. matrix_b_type == dbcsr_type_antisymmetric IF (matrix_a_symm .AND. matrix_b_symm) sym_fac = REAL(2.0, real_4) ! tracing a symmetric with a general matrix is not implemented, as it would require communication of blocks IF (matrix_a_symm .NEQV. matrix_b_symm) & DBCSR_ABORT("Tracing general with symmetric matrix NYI") a_row_blk_size => array_data(matrix_a%row_blk_size) a_col_blk_size => array_data(matrix_a%col_blk_size) b_row_blk_size => array_data(matrix_b%row_blk_size) b_col_blk_size => array_data(matrix_b%col_blk_size) CALL dbcsr_get_data(matrix_a%data_area, a_data) CALL dbcsr_get_data(matrix_b%data_area, b_data) ! let's go trace = REAL(0.0, real_4) IF (matrix_a%nblkrows_total .NE. matrix_b%nblkrows_total) & DBCSR_ABORT("this combination of transpose is NYI") DO row = 1, matrix_a%nblkrows_total a_row_size = a_row_blk_size(row) b_row_size = b_row_blk_size(row) IF (a_row_size .NE. b_row_size) DBCSR_ABORT("matrices not consistent") b_blk = matrix_b%row_p(row) + 1 b_frst_blk = matrix_b%row_p(row) + 1 b_last_blk = matrix_b%row_p(row + 1) DO a_blk = matrix_a%row_p(row) + 1, matrix_a%row_p(row + 1) IF (matrix_a%blk_p(a_blk) .EQ. 0) CYCLE ! Deleted block a_col = matrix_a%col_i(a_blk) a_col_size = a_col_blk_size(a_col) ! ! find the b_blk we assume here that the columns are ordered ! CALL dbcsr_find_column(a_col, b_frst_blk, b_last_blk, matrix_b%col_i, & matrix_b%blk_p, b_blk, found) IF (found) THEN b_col_size = b_col_blk_size(a_col) IF (a_col_size .NE. b_col_size) DBCSR_ABORT("matrices not consistent") ! nze = a_row_size*a_col_size ! IF (nze .GT. 0) THEN ! ! let's trace the blocks a_beg = ABS(matrix_a%blk_p(a_blk)) a_end = a_beg + nze - 1 b_beg = ABS(matrix_b%blk_p(b_blk)) b_end = b_beg + nze - 1 fac = REAL(1.0, real_4) IF (row .NE. a_col) fac = sym_fac trace = trace + fac*SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) END IF END IF END DO ! a_col END DO ! a_row ! ! sum CALL mp_sum(trace, dbcsr_mp_group(dbcsr_distribution_mp(matrix_a%dist))) END SUBROUTINE dbcsr_dot_c SUBROUTINE dbcsr_scale_c (matrix_a, alpha_scalar, last_column) !! Interface for matrix scaling by a scalar TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a COMPLEX(kind=real_4), INTENT(IN) :: alpha_scalar INTEGER, INTENT(IN), OPTIONAL :: last_column CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_scale_c' INTEGER :: error_handler TYPE(dbcsr_scalar_type) :: sc sc = dbcsr_scalar(alpha_scalar) CALL dbcsr_scalar_fill_all(sc) sc%data_type = dbcsr_get_data_type(matrix_a) CALL timeset(routineN, error_handler) IF (PRESENT(last_column)) THEN CALL dbcsr_scale_anytype(matrix_a, & alpha_scalar=sc, & limits=(/0, 0, 0, last_column/)) ELSE CALL dbcsr_scale_anytype(matrix_a, alpha_scalar=sc) END IF CALL timestop(error_handler) END SUBROUTINE dbcsr_scale_c SUBROUTINE dbcsr_scale_by_vector_c (matrix_a, alpha, side) !! Interface for matrix scaling by a vector TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a COMPLEX(kind=real_4), DIMENSION(:), INTENT(IN), TARGET, CONTIGUOUS :: alpha CHARACTER(LEN=*), INTENT(IN) :: side COMPLEX(kind=real_4), DIMENSION(:), POINTER, CONTIGUOUS :: tmp_p TYPE(dbcsr_data_obj) :: enc_alpha_vec CALL dbcsr_data_init(enc_alpha_vec) CALL dbcsr_data_new(enc_alpha_vec, dbcsr_type_complex_4) tmp_p => alpha CALL dbcsr_data_set_pointer(enc_alpha_vec, tmp_p) CALL dbcsr_scale_by_vector_anytype(matrix_a, enc_alpha_vec, side) CALL dbcsr_data_clear_pointer(enc_alpha_vec) CALL dbcsr_data_release(enc_alpha_vec) END SUBROUTINE dbcsr_scale_by_vector_c SUBROUTINE dbcsr_set_c (matrix, alpha) !! Interface for dbcsr_set TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_4), INTENT(IN) :: alpha CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set' INTEGER :: col, handle, row TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_4), DIMENSION(:, :), POINTER :: block LOGICAL :: tr CALL timeset(routineN, handle) IF (alpha == CMPLX(0.0, 0.0, real_4)) THEN CALL dbcsr_zero(matrix) ELSE IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_4) & DBCSR_ABORT("Incompatible data types") !TODO: could be speedup by direct assignment to data_area, similar to dbcsr_zero() CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, block, tr) block(:, :) = alpha END DO CALL dbcsr_iterator_stop(iter) END IF CALL timestop(handle) END SUBROUTINE dbcsr_set_c SUBROUTINE dbcsr_filter_c (matrix, eps, method, use_absolute, & filter_diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_4), INTENT(IN) :: eps INTEGER, INTENT(IN), OPTIONAL :: method LOGICAL, INTENT(in), OPTIONAL :: use_absolute, filter_diag CALL dbcsr_filter_anytype(matrix, dbcsr_scalar(eps), method, & use_absolute, filter_diag) END SUBROUTINE dbcsr_filter_c SUBROUTINE dbcsr_set_diag_c (matrix, diag) TYPE(dbcsr_type), INTENT(INOUT) :: matrix COMPLEX(kind=real_4), DIMENSION(:), INTENT(IN) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_set_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_4) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) block(i, i) = diag(row_offset + i - 1) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_set_diag_c SUBROUTINE dbcsr_get_diag_c (matrix, diag) TYPE(dbcsr_type), INTENT(IN) :: matrix COMPLEX(kind=real_4), DIMENSION(:), INTENT(OUT) :: diag CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_get_diag' INTEGER :: icol, irow, row_offset, handle, i LOGICAL :: tr TYPE(dbcsr_iterator) :: iter COMPLEX(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_4) & DBCSR_ABORT("Incompatible data types") IF (dbcsr_nfullrows_total(matrix) /= SIZE(diag)) & DBCSR_ABORT("Diagonal has wrong size") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") diag(:) = CMPLX(0.0, 0.0, real_4) CALL dbcsr_iterator_start(iter, matrix) DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, irow, icol, block, tr, row_offset=row_offset) IF (irow /= icol) CYCLE IF (sIZE(block, 1) /= sIZE(block, 2)) & DBCSR_ABORT("Diagonal block non-squared") DO i = 1, sIZE(block, 1) diag(row_offset + i - 1) = block(i, i) END DO END DO CALL dbcsr_iterator_stop(iter) CALL timestop(handle) END SUBROUTINE dbcsr_get_diag_c SUBROUTINE dbcsr_add_on_diag_c (matrix, alpha) !! add a constant to the diagonal of a matrix TYPE(dbcsr_type), INTENT(INOUT) :: matrix !! DBCSR matrix COMPLEX(kind=real_4), INTENT(IN) :: alpha !! scalar CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_add_on_diag' INTEGER :: handle, mynode, node, irow, i, row_size LOGICAL :: found, tr COMPLEX(kind=real_4), DIMENSION(:, :), POINTER :: block CALL timeset(routineN, handle) IF (dbcsr_get_data_type(matrix) /= dbcsr_type_complex_4) & DBCSR_ABORT("Incompatible data types") IF (.NOT. array_equality(matrix%row_blk_offset, matrix%col_blk_offset)) & DBCSR_ABORT("matrix not quadratic") mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix))) CALL dbcsr_work_create(matrix, work_mutable=.TRUE.) DO irow = 1, dbcsr_nblkrows_total(matrix) CALL dbcsr_get_stored_coordinates(matrix, irow, irow, node) IF (node /= mynode) CYCLE CALL dbcsr_get_block_p(matrix, irow, irow, block, tr, found, row_size=row_size) IF (.NOT. found) THEN ALLOCATE (block(row_size, row_size)) block(:, :) = CMPLX(0.0, 0.0, real_4) END IF DO i = 1, row_size block(i, i) = block(i, i) + alpha END DO IF (.NOT. found) THEN CALL dbcsr_put_block(matrix, irow, irow, block) DEALLOCATE (block) END IF END DO CALL dbcsr_finalize(matrix) CALL timestop(handle) END SUBROUTINE dbcsr_add_on_diag_c SUBROUTINE dbcsr_update_contiguous_blocks_c (matrix_a, matrix_b, first_lb_a, first_lb_b, nze, & do_scale, my_beta_scalar, found, iw) !! Low level function to sum contiguous chunks of blocks of the matrices (matrix_a = matrix_a + beta*matrix_b) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER, INTENT(IN) :: first_lb_a, first_lb_b, nze, iw LOGICAL, INTENT(IN) :: found, do_scale INTEGER :: ub_a, ub_b ub_a = first_lb_a + nze - 1 ub_b = first_lb_b + nze - 1 IF (found) THEN IF (do_scale) THEN CALL caxpy(nze, my_beta_scalar%c_sp, & matrix_b%data_area%d%c_sp (first_lb_b:ub_b), 1, & matrix_a%data_area%d%c_sp (first_lb_a:ub_a), 1) ELSE matrix_a%data_area%d%c_sp (first_lb_a:ub_a) = & matrix_a%data_area%d%c_sp (first_lb_a:ub_a) + & matrix_b%data_area%d%c_sp (first_lb_b:ub_b) END IF ELSE IF (do_scale) THEN matrix_a%wms(iw)%data_area%d%c_sp (first_lb_a:ub_a) = & my_beta_scalar%c_sp* & matrix_b%data_area%d%c_sp (first_lb_b:ub_b) ELSE matrix_a%wms(iw)%data_area%d%c_sp (first_lb_a:ub_a) = & matrix_b%data_area%d%c_sp (first_lb_b:ub_b) END IF END IF END SUBROUTINE dbcsr_update_contiguous_blocks_c SUBROUTINE dbcsr_add_anytype_c (matrix_a, matrix_b, iter, iw, do_scale, & my_beta_scalar, my_flop) !! Low level function to sum two matrices (matrix_a = matrix_a + beta*matrix_b TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a !! DBCSR matrix TYPE(dbcsr_type), INTENT(IN) :: matrix_b !! DBCSR matrix TYPE(dbcsr_iterator), INTENT(INOUT) :: iter INTEGER, INTENT(IN) :: iw LOGICAL, INTENT(IN) :: do_scale TYPE(dbcsr_scalar_type), INTENT(IN) :: my_beta_scalar INTEGER(KIND=int_8), INTENT(INOUT) :: my_flop INTEGER :: row, col, row_size, col_size, & nze, tot_nze, blk, & lb_a, first_lb_a, lb_a_val, & lb_b, first_lb_b INTEGER, DIMENSION(2) :: lb_row_blk LOGICAL :: was_found, found, tr ! some start values lb_row_blk(:) = 0 first_lb_a = matrix_a%wms(iw)%datasize + 1 first_lb_b = 0 tot_nze = 0 ! DO WHILE (dbcsr_iterator_blocks_left(iter)) CALL dbcsr_iterator_next_block(iter, row, col, blk, tr, lb_b, row_size, col_size) nze = row_size*col_size IF (nze .LE. 0) CYCLE IF (lb_row_blk(1) .LT. row) THEN lb_row_blk(1) = row lb_row_blk(2) = matrix_a%row_p(row) + 1 END IF ! get b-block index lb_b = ABS(lb_b) CALL dbcsr_find_column(col, lb_row_blk(2), matrix_a%row_p(row + 1), matrix_a%col_i, matrix_a%blk_p, blk, found) lb_row_blk(2) = blk + 1 ! get index of a-block lb_a whether found (from matrix_a) or not (from workspace array) IF (found) THEN my_flop = my_flop + nze*2 lb_a = ABS(matrix_a%blk_p(blk)) ELSE lb_a = matrix_a%wms(iw)%datasize + 1 lb_a_val = lb_a IF (tr) lb_a_val = -lb_a matrix_a%wms(iw)%lastblk = matrix_a%wms(iw)%lastblk + 1 matrix_a%wms(iw)%row_i(matrix_a%wms(iw)%lastblk) = row matrix_a%wms(iw)%col_i(matrix_a%wms(iw)%lastblk) = col matrix_a%wms(iw)%blk_p(matrix_a%wms(iw)%lastblk) = lb_a_val matrix_a%wms(iw)%datasize = matrix_a%wms(iw)%datasize + nze END IF ! at the first iteration we skip this and go directly to initialization after IF (first_lb_b .NE. 0) THEN ! if found status is the same as before then probably we are in contiguous blocks IF ((found .EQV. was_found) .AND. & (first_lb_b + tot_nze .EQ. lb_b) .AND. & (first_lb_a + tot_nze) .EQ. lb_a) THEN tot_nze = tot_nze + nze CYCLE END IF ! save block chunk CALL dbcsr_update_contiguous_blocks_c (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF ! first_lb_a = lb_a first_lb_b = lb_b tot_nze = nze was_found = found END DO ! save the last block or chunk of blocks IF (first_lb_b .NE. 0) THEN call dbcsr_update_contiguous_blocks_c (matrix_a, matrix_b, first_lb_a, first_lb_b, tot_nze, & do_scale, my_beta_scalar, was_found, iw) END IF END SUBROUTINE dbcsr_add_anytype_c # 3132 "/__w/dbcsr/dbcsr/src/ops/dbcsr_operations.F" END MODULE dbcsr_operations