diff --git a/src/core/dbcsr_types.F b/src/core/dbcsr_types.F index a2e17ad8b05..1f01367934a 100644 --- a/src/core/dbcsr_types.F +++ b/src/core/dbcsr_types.F @@ -84,6 +84,7 @@ MODULE dbcsr_types PUBLIC :: dbcsr_type_invalid, dbcsr_type_no_symmetry, dbcsr_type_symmetric, & dbcsr_type_antisymmetric, dbcsr_type_hermitian, dbcsr_type_antihermitian + PUBLIC :: dbcsr_symmetry_data_type_compatible PUBLIC :: dbcsr_no_transpose, dbcsr_transpose, dbcsr_conjugate_transpose PUBLIC :: dbcsr_repl_none, dbcsr_repl_row, dbcsr_repl_col, dbcsr_repl_full @@ -588,4 +589,28 @@ MODULE dbcsr_types INTEGER(KIND=int_8), DIMENSION(SIZE(dbcsr_mpi_size_limits) + 1, 2, 2) :: data_size_breakdown = -1 END TYPE dbcsr_mpi_statistics_type +CONTAINS + + PURE FUNCTION dbcsr_symmetry_data_type_compatible(matrix_type, data_type) RESULT(compatible) + !! Whether a matrix symmetry is supported for the given scalar type. + + CHARACTER, INTENT(IN) :: matrix_type + INTEGER, INTENT(IN) :: data_type + LOGICAL :: compatible + + compatible = .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) + compatible = .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) + compatible = .TRUE. + END SELECT + END SELECT + END FUNCTION dbcsr_symmetry_data_type_compatible + END MODULE dbcsr_types diff --git a/src/dbcsr_api.F b/src/dbcsr_api.F index a0a3faa4d10..4e20d77e8ee 100644 --- a/src/dbcsr_api.F +++ b/src/dbcsr_api.F @@ -123,9 +123,10 @@ MODULE dbcsr_api dbcsr_func_inverse, dbcsr_func_tanh, dbcsr_iterator_prv => dbcsr_iterator, dbcsr_mp_obj, & dbcsr_no_transpose, dbcsr_norm_column, dbcsr_norm_frobenius, dbcsr_norm_maxabsnorm, & dbcsr_prv_type => dbcsr_type, dbcsr_scalar_type, dbcsr_transpose, & - dbcsr_type_antisymmetric, dbcsr_type_complex_4, dbcsr_type_complex_8, & + dbcsr_conjugate_transpose, dbcsr_type_antihermitian, dbcsr_type_antisymmetric, & + dbcsr_type_complex_4, dbcsr_type_complex_8, & dbcsr_type_complex_default, dbcsr_type_no_symmetry, dbcsr_type_real_4, dbcsr_type_real_8, & - dbcsr_type_real_default, dbcsr_type_symmetric + dbcsr_type_real_default, dbcsr_type_symmetric, dbcsr_type_hermitian USE dbcsr_dist_util, ONLY: dbcsr_convert_offsets_to_sizes => convert_offsets_to_sizes, & dbcsr_convert_sizes_to_offsets => convert_sizes_to_offsets, & dbcsr_checksum_prv => dbcsr_checksum, & @@ -150,7 +151,10 @@ MODULE dbcsr_api PUBLIC :: dbcsr_type_no_symmetry PUBLIC :: dbcsr_type_symmetric PUBLIC :: dbcsr_type_antisymmetric + PUBLIC :: dbcsr_type_hermitian + PUBLIC :: dbcsr_type_antihermitian PUBLIC :: dbcsr_transpose + PUBLIC :: dbcsr_conjugate_transpose PUBLIC :: dbcsr_no_transpose PUBLIC :: dbcsr_type_complex_8 PUBLIC :: dbcsr_type_real_4 @@ -1401,11 +1405,14 @@ SUBROUTINE dbcsr_trace_${nametype1}$ (matrix_a, trace) CALL dbcsr_trace_prv(matrix_a%prv, trace) END SUBROUTINE dbcsr_trace_${nametype1}$ - SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, result) + SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, result, conjugate_a) + !! Computes sum(A(i,j)*B(i,j)) by default. If conjugate_a is true, + !! computes the Frobenius inner product sum(conjg(A(i,j))*B(i,j)). TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b ${type1}$, INTENT(INOUT) :: result + LOGICAL, INTENT(IN), OPTIONAL :: conjugate_a - CALL dbcsr_dot_prv(matrix_a%prv, matrix_b%prv, result) + CALL dbcsr_dot_prv(matrix_a%prv, matrix_b%prv, result, conjugate_a) END SUBROUTINE dbcsr_dot_${nametype1}$ SUBROUTINE dbcsr_multiply_${nametype1}$ (transa, transb, & diff --git a/src/ops/dbcsr_operations.F b/src/ops/dbcsr_operations.F index 1c81facb86e..6fce739ad3c 100644 --- a/src/ops/dbcsr_operations.F +++ b/src/ops/dbcsr_operations.F @@ -77,7 +77,7 @@ MODULE dbcsr_operations 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 + dbcsr_type_real_8, dbcsr_type_symmetric, dbcsr_symmetry_data_type_compatible USE dbcsr_dist_util, ONLY: find_block_of_element USE dbcsr_work_operations, ONLY: dbcsr_create, & dbcsr_finalize, & @@ -612,12 +612,17 @@ SUBROUTINE dbcsr_add_anytype(matrix_a, matrix_b, alpha_scalar, beta_scalar, flop 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 + IF (dbcsr_has_symmetry(matrix_a) .NEQV. dbcsr_has_symmetry(matrix_b)) THEN DBCSR_ABORT("Summing general with symmetric matrix NYI") END IF + IF (dbcsr_has_symmetry(matrix_a) .AND. & + dbcsr_get_matrix_type(matrix_a) .NE. dbcsr_get_matrix_type(matrix_b)) & + DBCSR_ABORT("Summing matrices with different symmetry types NYI") + IF (.NOT. dbcsr_symmetry_data_type_compatible(dbcsr_get_matrix_type(matrix_a), & + dbcsr_get_data_type(matrix_a)) .OR. & + .NOT. dbcsr_symmetry_data_type_compatible(dbcsr_get_matrix_type(matrix_b), & + dbcsr_get_data_type(matrix_b))) & + DBCSR_ABORT("Matrix symmetry is incompatible with its scalar data type.") data_type_a = dbcsr_get_data_type(matrix_a) data_type_b = dbcsr_get_data_type(matrix_b) @@ -1180,22 +1185,7 @@ LOGICAL FUNCTION symmetry_consistent(matrix_type, data_type) 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 + symmetry_consistent = dbcsr_symmetry_data_type_compatible(matrix_type, data_type) END FUNCTION symmetry_consistent @@ -1246,7 +1236,7 @@ SUBROUTINE dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, & !! 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 + !! matrix type compatible with the scalar data type: N/S/A for real, N/H/K for complex CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_copy' CHARACTER :: new_matrix_type, repl_type @@ -1358,6 +1348,14 @@ SUBROUTINE dbcsr_copy_into_existing(matrix_b, matrix_a) 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. symmetry_consistent(dbcsr_get_matrix_type(matrix_a), data_type) .OR. & + .NOT. symmetry_consistent(dbcsr_get_matrix_type(matrix_b), data_type)) & + DBCSR_ABORT("Matrix symmetry is incompatible with its scalar data type.") + IF (dbcsr_has_symmetry(matrix_a) .AND. dbcsr_has_symmetry(matrix_b) .AND. & + dbcsr_get_matrix_type(matrix_a) .NE. dbcsr_get_matrix_type(matrix_b)) & + DBCSR_ABORT("Copying between different symmetric matrix types NYI.") + IF (.NOT. dbcsr_has_symmetry(matrix_b) .AND. dbcsr_has_symmetry(matrix_a)) & + DBCSR_ABORT("Copying a compressed matrix into an uncompressed target NYI.") neg_real = matrix_b%negate_real neg_imag = matrix_b%negate_imaginary making_symmetric = dbcsr_has_symmetry(matrix_b) & @@ -2189,8 +2187,7 @@ FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm) 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 + any_sym = dbcsr_has_symmetry(matrix) IF (nr .NE. nc) & DBCSR_ABORT("not a square matrix") @@ -2228,8 +2225,7 @@ FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm) 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)) + 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) @@ -2239,8 +2235,7 @@ FUNCTION dbcsr_gershgorin_norm(matrix) RESULT(norm) 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)) + buff_d(col_offset + j - 1) = buff_d(col_offset + j - 1) + ABS(data_z(i, j)) END DO END DO CASE DEFAULT @@ -2325,8 +2320,7 @@ FUNCTION dbcsr_frobenius_norm(matrix, local) RESULT(norm) 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 + any_sym = dbcsr_has_symmetry(matrix) norm = 0.0_dp CALL dbcsr_iterator_start(iter, matrix) @@ -2345,14 +2339,12 @@ FUNCTION dbcsr_frobenius_norm(matrix, local) RESULT(norm) 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") + IF (any_sym .AND. row .NE. col) fac = 2.0_dp 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") + IF (any_sym .AND. row .NE. col) fac = 2.0_dp norm = norm + fac*REAL(SUM(CONJG(data_z)*data_z), KIND=real_8) CASE DEFAULT DBCSR_ABORT("Wrong data type") @@ -2591,7 +2583,7 @@ SUBROUTINE dbcsr_trace_sd(matrix_a, trace) CALL timestop(handle) END SUBROUTINE dbcsr_trace_sd - SUBROUTINE dbcsr_dot_sd(matrix_a, matrix_b, trace) + SUBROUTINE dbcsr_dot_sd(matrix_a, matrix_b, trace, conjugate_a) !! Dot product of DBCSR matrices !! \result the dot product of the matrices @@ -2599,6 +2591,7 @@ SUBROUTINE dbcsr_dot_sd(matrix_a, matrix_b, trace) !! DBCSR matrices !! DBCSR matrices REAL(kind=real_8), INTENT(INOUT) :: trace + LOGICAL, INTENT(IN), OPTIONAL :: conjugate_a CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_dot_sd' @@ -2608,11 +2601,11 @@ SUBROUTINE dbcsr_dot_sd(matrix_a, matrix_b, trace) 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) + CALL dbcsr_dot_d(matrix_a, matrix_b, trace, conjugate_a) 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) + CALL dbcsr_dot_s(matrix_a, matrix_b, trace_4, conjugate_a) trace = REAL(trace_4, real_8) ELSE DBCSR_ABORT("Invalid combination of data type, NYI") @@ -2689,14 +2682,17 @@ SUBROUTINE dbcsr_trace_${nametype1}$ (matrix_a, trace) CALL timestop(error_handle) END SUBROUTINE dbcsr_trace_${nametype1}$ - SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, trace) + SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, trace, conjugate_a) !! Dot product of DBCSR matrices + !! By default this is the bilinear product sum(A(i,j)*B(i,j)). Set + !! conjugate_a to compute the Frobenius inner product sum(conjg(A(i,j))*B(i,j)). TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b !! DBCSR matrices !! DBCSR matrices ${type1}$, INTENT(INOUT) :: trace !! the trace of the product of the matrices + LOGICAL, INTENT(IN), OPTIONAL :: conjugate_a 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 @@ -2704,8 +2700,8 @@ SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, trace) INTEGER, DIMENSION(:), POINTER :: a_col_blk_size, & a_row_blk_size, & b_col_blk_size, b_row_blk_size - ${type1}$ :: sym_fac, fac - LOGICAL :: found, matrix_a_symm, matrix_b_symm + ${type1}$ :: dot_block, sym_fac + LOGICAL :: found, matrix_a_symm, matrix_b_symm, my_conjugate_a ${type1}$, DIMENSION(:), POINTER :: a_data, b_data ! --------------------------------------------------------------------------- @@ -2714,13 +2710,15 @@ SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, trace) .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, ${kind1}$) 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, ${kind1}$) + matrix_a_symm = matrix_a_type /= dbcsr_type_no_symmetry + matrix_b_symm = matrix_b_type /= dbcsr_type_no_symmetry + sym_fac = REAL(1.0, ${kind1}$) + IF (matrix_a_symm .AND. matrix_b_symm .AND. matrix_a_type /= matrix_b_type) & + sym_fac = REAL(-1.0, ${kind1}$) + my_conjugate_a = .FALSE. + IF (PRESENT(conjugate_a)) my_conjugate_a = conjugate_a ! 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) & @@ -2766,10 +2764,25 @@ SUBROUTINE dbcsr_dot_${nametype1}$ (matrix_a, matrix_b, trace) 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, ${kind1}$) - IF (row .NE. a_col) fac = sym_fac - - trace = trace + fac*SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) + #:if n >= 2 + IF (my_conjugate_a) THEN + dot_block = SUM(CONJG(a_data(a_beg:a_end))*b_data(b_beg:b_end)) + ELSE + dot_block = SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) + END IF + IF (matrix_a_symm .AND. row .NE. a_col) THEN + trace = trace + dot_block + sym_fac*CONJG(dot_block) + ELSE + trace = trace + dot_block + END IF + #:else + dot_block = SUM(a_data(a_beg:a_end)*b_data(b_beg:b_end)) + IF (matrix_a_symm .AND. row .NE. a_col) THEN + trace = trace + (1.0_${kind1}$+sym_fac)*dot_block + ELSE + trace = trace + dot_block + END IF + #:endif END IF END IF diff --git a/src/work/dbcsr_work_operations.F b/src/work/dbcsr_work_operations.F index 5a441927b77..e4eebdc4426 100644 --- a/src/work/dbcsr_work_operations.F +++ b/src/work/dbcsr_work_operations.F @@ -78,7 +78,8 @@ MODULE dbcsr_work_operations dbcsr_slot_nfullrows_total, dbcsr_slot_nze, dbcsr_slot_row_p, dbcsr_slot_size, 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_real_default, dbcsr_type_symmetric, dbcsr_work_type + dbcsr_type_real_8, dbcsr_type_real_default, dbcsr_type_symmetric, & + dbcsr_symmetry_data_type_compatible, dbcsr_work_type USE dbcsr_dist_util, ONLY: convert_sizes_to_offsets, & dbcsr_calc_block_sizes, & dbcsr_verify_matrix, & @@ -134,7 +135,8 @@ SUBROUTINE dbcsr_create_new(matrix, name, dist, matrix_type, & TYPE(dbcsr_distribution_obj), INTENT(IN) :: dist !! distribution_2d distribution CHARACTER, INTENT(IN) :: matrix_type - !! 'N' for normal, 'T' for transposed, 'S' for symmetric, and 'A' for antisymmetric + !! Matrix storage symmetry: 'N' is uncompressed. Real data supports 'S' and 'A'; + !! complex data supports 'H' and 'K'. Other scalar/symmetry combinations are rejected. INTEGER, DIMENSION(:), INTENT(INOUT), POINTER, & CONTIGUOUS, OPTIONAL :: row_blk_size, col_blk_size TYPE(array_i1d_obj), INTENT(IN), OPTIONAL :: row_blk_size_obj, col_blk_size_obj @@ -158,7 +160,7 @@ SUBROUTINE dbcsr_create_new(matrix, name, dist, matrix_type, & CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_create_new' CHARACTER :: matrix_type_l - INTEGER :: handle, my_nze + INTEGER :: handle, matrix_data_type, my_nze INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: vec_col_blk_offset, vec_row_blk_offset INTEGER, DIMENSION(dbcsr_meta_size) :: new_meta LOGICAL :: hijack, my_make_index @@ -169,6 +171,14 @@ SUBROUTINE dbcsr_create_new(matrix, name, dist, matrix_type, & CALL timeset(routineN, handle) + matrix_type_l = matrix_type + CALL uppercase(matrix_type_l) + matrix_data_type = dbcsr_type_real_default + IF (PRESENT(data_type)) matrix_data_type = data_type + IF (.NOT. dbcsr_symmetry_data_type_compatible(matrix_type_l, matrix_data_type)) & + DBCSR_ABORT("Matrix symmetry '"//matrix_type_l// & + "' is incompatible with its scalar data type.") + ! Reuse matrix only if has actually been allocated. hijack = ASSOCIATED(matrix%index) IF (PRESENT(reuse)) hijack = reuse @@ -397,8 +407,6 @@ SUBROUTINE dbcsr_create_new(matrix, name, dist, matrix_type, & matrix%negate_real = .FALSE. matrix%negate_imaginary = .FALSE. !matrix%transpose = .FALSE. - matrix_type_l = matrix_type - CALL uppercase(matrix_type_l) SELECT CASE (matrix_type_l) CASE (dbcsr_type_no_symmetry) CASE (dbcsr_type_symmetric) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 172dcd7146a..40a88127d45 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -85,7 +85,7 @@ if (NOT (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray")) endif () # For each test, set a variable testname_SRCS defining the sources of that test -set(dbcsr_unittest1_SRCS dbcsr_unittest1.F) +set(dbcsr_unittest1_SRCS dbcsr_unittest1.F dbcsr_test_symmetry.F) set(dbcsr_unittest2_SRCS dbcsr_unittest2.F) set(dbcsr_unittest3_SRCS dbcsr_unittest3.F) set(dbcsr_unittest4_SRCS dbcsr_unittest4.F dbcsr_test_scale_by_vector.F) diff --git a/tests/dbcsr_test_add.F b/tests/dbcsr_test_add.F index 7366a76c2ae..981f87a4d88 100644 --- a/tests/dbcsr_test_add.F +++ b/tests/dbcsr_test_add.F @@ -41,7 +41,8 @@ MODULE dbcsr_test_add dbcsr_type_antihermitian, dbcsr_type_antisymmetric, dbcsr_type_complex_4, & dbcsr_type_complex_4_2d, dbcsr_type_complex_8, dbcsr_type_complex_8_2d, & dbcsr_type_hermitian, dbcsr_type_no_symmetry, dbcsr_type_real_4, dbcsr_type_real_4_2d, & - dbcsr_type_real_8, dbcsr_type_real_8_2d, dbcsr_type_symmetric + dbcsr_type_real_8, dbcsr_type_real_8_2d, dbcsr_type_symmetric, & + dbcsr_symmetry_data_type_compatible USE dbcsr_work_operations, ONLY: dbcsr_create #include "base/dbcsr_base_uses.f90" @@ -142,6 +143,8 @@ SUBROUTINE dbcsr_test_adds(test_name, mp_group, mp_env, npdims, io_unit, & TYPE = types(itype) do_complex = TYPE .EQ. dbcsr_type_complex_4 .OR. TYPE .EQ. dbcsr_type_complex_8 + IF (.NOT. dbcsr_symmetry_data_type_compatible(a_symm, TYPE) .OR. & + .NOT. dbcsr_symmetry_data_type_compatible(b_symm, TYPE)) CYCLE alpha_obj = compx_to_dbcsr_scalar(alpha, TYPE) beta_obj = compx_to_dbcsr_scalar(beta, TYPE) diff --git a/tests/dbcsr_test_multiply.F b/tests/dbcsr_test_multiply.F index 10e2b540c29..5602119f0b2 100644 --- a/tests/dbcsr_test_multiply.F +++ b/tests/dbcsr_test_multiply.F @@ -45,7 +45,8 @@ MODULE dbcsr_test_multiply dbcsr_no_transpose, dbcsr_scalar_type, dbcsr_transpose, dbcsr_type, & dbcsr_type_antisymmetric, dbcsr_type_complex_4, dbcsr_type_complex_4_2d, & dbcsr_type_complex_8, dbcsr_type_complex_8_2d, dbcsr_type_no_symmetry, dbcsr_type_real_4, & - dbcsr_type_real_4_2d, dbcsr_type_real_8, dbcsr_type_real_8_2d, dbcsr_type_symmetric + dbcsr_type_real_4_2d, dbcsr_type_real_8, dbcsr_type_real_8_2d, dbcsr_type_symmetric, & + dbcsr_symmetry_data_type_compatible USE dbcsr_work_operations, ONLY: dbcsr_create #include "base/dbcsr_base_uses.f90" @@ -171,6 +172,9 @@ SUBROUTINE dbcsr_test_multiplies(test_name, mp_group, mp_env, npdims, io_unit, & TYPE = types(itype) do_complex = TYPE .EQ. dbcsr_type_complex_4 .OR. TYPE .EQ. dbcsr_type_complex_8 + IF (.NOT. dbcsr_symmetry_data_type_compatible(a_symm, TYPE) .OR. & + .NOT. dbcsr_symmetry_data_type_compatible(b_symm, TYPE) .OR. & + .NOT. dbcsr_symmetry_data_type_compatible(c_symm, TYPE)) CYCLE alpha_obj = compx_to_dbcsr_scalar(alpha, TYPE) beta_obj = compx_to_dbcsr_scalar(beta, TYPE) diff --git a/tests/dbcsr_test_symmetry.F b/tests/dbcsr_test_symmetry.F new file mode 100644 index 00000000000..4741f22516d --- /dev/null +++ b/tests/dbcsr_test_symmetry.F @@ -0,0 +1,262 @@ +!--------------------------------------------------------------------------------------------------! +! 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_test_symmetry + USE dbcsr_api, ONLY: dbcsr_binary_read, dbcsr_binary_write, dbcsr_conjugate_transpose, & + dbcsr_copy, dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, & + dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_finalize, & + dbcsr_frobenius_norm, dbcsr_gershgorin_norm, dbcsr_get_data_type, & + dbcsr_get_matrix_type, dbcsr_multiply, dbcsr_no_transpose, & + dbcsr_put_block, dbcsr_release, dbcsr_transposed, 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, dbcsr_work_create + USE dbcsr_kinds, ONLY: real_8 + USE dbcsr_mpiwrap, ONLY: mp_cart_rank, mp_comm_type + USE dbcsr_types, ONLY: dbcsr_symmetry_data_type_compatible + IMPLICIT NONE + PRIVATE + + PUBLIC :: dbcsr_test_hermitian_symmetry + +CONTAINS + + LOGICAL FUNCTION dbcsr_test_hermitian_symmetry(group, npdims, io_unit) RESULT(success) + TYPE(mp_comm_type), INTENT(IN) :: group + INTEGER, DIMENSION(2), INTENT(IN) :: npdims + INTEGER, INTENT(IN) :: io_unit + + INTEGER, DIMENSION(:), POINTER :: col_dist, row_dist, col_blk_size, row_blk_size + INTEGER, DIMENSION(:, :), POINTER :: pgrid + INTEGER :: pcol, prow + COMPLEX(real_8) :: dot_value + TYPE(dbcsr_distribution_type) :: dist + TYPE(dbcsr_type) :: h, k, identity, h_copy, k_copy, h_transposed, k_transposed, & + h_product, k_product, h_read, k_read + + success = .TRUE. + CALL check(dbcsr_api_type_exports(), "high-level symmetry and transpose exports", success, io_unit) + CALL check_type_combinations(success, io_unit) + + ALLOCATE (pgrid(0:npdims(1) - 1, 0:npdims(2) - 1)) + DO prow = 0, npdims(1) - 1 + DO pcol = 0, npdims(2) - 1 + CALL mp_cart_rank(group, (/prow, pcol/), pgrid(prow, pcol)) + END DO + END DO + ALLOCATE (row_dist(2), col_dist(2)) + row_dist = 0 + col_dist = 0 + CALL dbcsr_distribution_new(dist, group=group%get_handle(), pgrid=pgrid, & + row_dist=row_dist, col_dist=col_dist) + DEALLOCATE (pgrid, row_dist, col_dist) + + CALL make_complex_matrix(h, dist, dbcsr_type_hermitian, & + CMPLX(2.0_real_8, 0.0_real_8, KIND=real_8), & + CMPLX(1.0_real_8, 2.0_real_8, KIND=real_8), & + CMPLX(3.0_real_8, 0.0_real_8, KIND=real_8)) + CALL make_complex_matrix(k, dist, dbcsr_type_antihermitian, & + CMPLX(0.0_real_8, 2.0_real_8, KIND=real_8), & + CMPLX(1.0_real_8, 2.0_real_8, KIND=real_8), & + CMPLX(0.0_real_8, -3.0_real_8, KIND=real_8)) + CALL make_complex_matrix(identity, dist, dbcsr_type_no_symmetry, & + CMPLX(1.0_real_8, 0.0_real_8, KIND=real_8), & + CMPLX(0.0_real_8, 0.0_real_8, KIND=real_8), & + CMPLX(1.0_real_8, 0.0_real_8, KIND=real_8)) + + CALL check(dbcsr_get_matrix_type(h) == dbcsr_type_hermitian .AND. & + dbcsr_get_data_type(h) == dbcsr_type_complex_8, "Hermitian creation", success, io_unit) + CALL check(dbcsr_get_matrix_type(k) == dbcsr_type_antihermitian .AND. & + dbcsr_get_data_type(k) == dbcsr_type_complex_8, "anti-Hermitian creation", success, io_unit) + + CALL dbcsr_dot(h, h, dot_value) + CALL check_complex(dot_value, CMPLX(7.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian bilinear dot", success, io_unit) + CALL dbcsr_dot(h, h, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian Frobenius product", success, io_unit) + CALL dbcsr_dot(k, k, dot_value) + CALL check_complex(dot_value, CMPLX(-19.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian bilinear dot", success, io_unit) + CALL dbcsr_dot(k, k, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian Frobenius product", success, io_unit) + CALL dbcsr_dot(h, k, dot_value) + CALL check_complex(dot_value, CMPLX(0.0_real_8, 3.0_real_8, KIND=real_8), & + "mixed H/K bilinear dot", success, io_unit) + CALL dbcsr_dot(h, k, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(0.0_real_8, -5.0_real_8, KIND=real_8), & + "mixed H/K Frobenius product", success, io_unit) + + CALL check_real(dbcsr_frobenius_norm(h), SQRT(23.0_real_8), & + "Hermitian Frobenius norm", success, io_unit) + CALL check_real(dbcsr_frobenius_norm(k), SQRT(23.0_real_8), & + "anti-Hermitian Frobenius norm", success, io_unit) + CALL check_real(dbcsr_gershgorin_norm(h), 3.0_real_8 + SQRT(5.0_real_8), & + "Hermitian Gershgorin norm", success, io_unit) + CALL check_real(dbcsr_gershgorin_norm(k), 3.0_real_8 + SQRT(5.0_real_8), & + "anti-Hermitian Gershgorin norm", success, io_unit) + + CALL dbcsr_copy(h_copy, h) + CALL dbcsr_copy(k_copy, k) + CALL check(dbcsr_get_matrix_type(h_copy) == dbcsr_type_hermitian .AND. & + dbcsr_get_matrix_type(k_copy) == dbcsr_type_antihermitian, & + "copy preserves H/K symmetry", success, io_unit) + + CALL dbcsr_transposed(h_transposed, h, shallow_data_copy=.FALSE., transpose_distribution=.FALSE.) + CALL dbcsr_transposed(k_transposed, k, shallow_data_copy=.FALSE., transpose_distribution=.FALSE.) + CALL dbcsr_dot(h_transposed, h, dot_value) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian transpose", success, io_unit) + CALL dbcsr_dot(k_transposed, k, dot_value) + CALL check_complex(dot_value, CMPLX(-23.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian transpose", success, io_unit) + + CALL dbcsr_filter(h, 0.1_real_8) + CALL dbcsr_filter(k, 0.1_real_8) + CALL dbcsr_dot(h, h, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian filter preserves stored values", success, io_unit) + CALL dbcsr_dot(k, k, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian filter preserves stored values", success, io_unit) + + CALL dbcsr_multiply(dbcsr_no_transpose, dbcsr_no_transpose, & + CMPLX(1.0_real_8, 0.0_real_8, KIND=real_8), h, identity, & + CMPLX(0.0_real_8, 0.0_real_8, KIND=real_8), h_product) + CALL dbcsr_multiply(dbcsr_conjugate_transpose, dbcsr_no_transpose, & + CMPLX(-1.0_real_8, 0.0_real_8, KIND=real_8), k, identity, & + CMPLX(0.0_real_8, 0.0_real_8, KIND=real_8), k_product) + CALL check(dbcsr_get_matrix_type(h_product) == dbcsr_type_hermitian .AND. & + dbcsr_get_matrix_type(k_product) == dbcsr_type_antihermitian, & + "multiplication preserves requested H/K output types", success, io_unit) + CALL dbcsr_dot(h_product, h, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian multiplication and conjugate transpose", success, io_unit) + CALL dbcsr_dot(k_product, k, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian multiplication and conjugate transpose", success, io_unit) + + CALL dbcsr_binary_write(h, "dbcsr_test_hermitian_io.bin") + CALL dbcsr_binary_read("dbcsr_test_hermitian_io.bin", dist, h_read) + CALL check(dbcsr_get_matrix_type(h_read) == dbcsr_type_hermitian, & + "Hermitian binary read preserves symmetry", success, io_unit) + CALL dbcsr_dot(h_read, h, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "Hermitian binary round trip", success, io_unit) + CALL dbcsr_binary_write(k, "dbcsr_test_antihermitian_io.bin") + CALL dbcsr_binary_read("dbcsr_test_antihermitian_io.bin", dist, k_read) + CALL check(dbcsr_get_matrix_type(k_read) == dbcsr_type_antihermitian, & + "anti-Hermitian binary read preserves symmetry", success, io_unit) + CALL dbcsr_dot(k_read, k, dot_value, conjugate_a=.TRUE.) + CALL check_complex(dot_value, CMPLX(23.0_real_8, 0.0_real_8, KIND=real_8), & + "anti-Hermitian binary round trip", success, io_unit) + + CALL dbcsr_release(h_read) + CALL dbcsr_release(k_read) + CALL dbcsr_release(h_product) + CALL dbcsr_release(k_product) + CALL dbcsr_release(h_transposed) + CALL dbcsr_release(k_transposed) + CALL dbcsr_release(h_copy) + CALL dbcsr_release(k_copy) + CALL dbcsr_release(identity) + CALL dbcsr_release(h) + CALL dbcsr_release(k) + CALL dbcsr_distribution_release(dist) + END FUNCTION dbcsr_test_hermitian_symmetry + + SUBROUTINE make_complex_matrix(matrix, dist, matrix_type, diag1, offdiag, diag2) + TYPE(dbcsr_type), INTENT(INOUT) :: matrix + TYPE(dbcsr_distribution_type), INTENT(IN):: dist + CHARACTER, INTENT(IN) :: matrix_type + COMPLEX(real_8), INTENT(IN) :: diag1, offdiag, diag2 + + INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size + COMPLEX(real_8), DIMENSION(1, 1) :: block + + ALLOCATE (row_blk_size(2), col_blk_size(2)) + row_blk_size = 1 + col_blk_size = 1 + CALL dbcsr_create(matrix, "symmetry test", dist, matrix_type, row_blk_size, col_blk_size, & + data_type=dbcsr_type_complex_8) + DEALLOCATE (row_blk_size, col_blk_size) + CALL dbcsr_work_create(matrix, nblks_guess=3, sizedata_guess=3) + block(1, 1) = diag1 + CALL dbcsr_put_block(matrix, 1, 1, block) + block(1, 1) = offdiag + IF (ABS(offdiag) > 0.0_real_8) CALL dbcsr_put_block(matrix, 1, 2, block) + block(1, 1) = diag2 + CALL dbcsr_put_block(matrix, 2, 2, block) + CALL dbcsr_finalize(matrix) + END SUBROUTINE make_complex_matrix + + PURE LOGICAL FUNCTION dbcsr_api_type_exports() RESULT(ok) + USE dbcsr_api, ONLY: api_conjugate_transpose => dbcsr_conjugate_transpose, & + api_antihermitian => dbcsr_type_antihermitian, & + api_hermitian => dbcsr_type_hermitian + ok = api_hermitian == 'H' .AND. api_antihermitian == 'K' .AND. api_conjugate_transpose == 'C' + END FUNCTION dbcsr_api_type_exports + + SUBROUTINE check_type_combinations(success, io_unit) + LOGICAL, INTENT(INOUT) :: success + INTEGER, INTENT(IN) :: io_unit + + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_no_symmetry, dbcsr_type_real_4), & + "real single N combination", success, io_unit) + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_symmetric, dbcsr_type_real_8), & + "real double S combination", success, io_unit) + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_antisymmetric, dbcsr_type_real_8), & + "real double A combination", success, io_unit) + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_no_symmetry, dbcsr_type_complex_4), & + "complex single N combination", success, io_unit) + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_hermitian, dbcsr_type_complex_8), & + "complex double H combination", success, io_unit) + CALL check(dbcsr_symmetry_data_type_compatible(dbcsr_type_antihermitian, dbcsr_type_complex_8), & + "complex double K combination", success, io_unit) + CALL check(.NOT. dbcsr_symmetry_data_type_compatible(dbcsr_type_hermitian, dbcsr_type_real_8), & + "real H combination rejected", success, io_unit) + CALL check(.NOT. dbcsr_symmetry_data_type_compatible(dbcsr_type_antihermitian, dbcsr_type_real_8), & + "real K combination rejected", success, io_unit) + CALL check(.NOT. dbcsr_symmetry_data_type_compatible(dbcsr_type_symmetric, dbcsr_type_complex_8), & + "complex S combination rejected", success, io_unit) + CALL check(.NOT. dbcsr_symmetry_data_type_compatible(dbcsr_type_antisymmetric, dbcsr_type_complex_8), & + "complex A combination rejected", success, io_unit) + END SUBROUTINE check_type_combinations + + SUBROUTINE check(condition, label, success, io_unit) + LOGICAL, INTENT(IN) :: condition + CHARACTER(LEN=*), INTENT(IN) :: label + LOGICAL, INTENT(INOUT) :: success + INTEGER, INTENT(IN) :: io_unit + + IF (condition) RETURN + success = .FALSE. + IF (io_unit > 0) WRITE (io_unit, *) "FAILED: ", TRIM(label) + END SUBROUTINE check + + SUBROUTINE check_complex(actual, expected, label, success, io_unit) + COMPLEX(real_8), INTENT(IN) :: actual, expected + CHARACTER(LEN=*), INTENT(IN) :: label + LOGICAL, INTENT(INOUT) :: success + INTEGER, INTENT(IN) :: io_unit + + CALL check(ABS(actual - expected) < 1.0E-12_real_8, label, success, io_unit) + END SUBROUTINE check_complex + + SUBROUTINE check_real(actual, expected, label, success, io_unit) + REAL(real_8), INTENT(IN) :: actual, expected + CHARACTER(LEN=*), INTENT(IN) :: label + LOGICAL, INTENT(INOUT) :: success + INTEGER, INTENT(IN) :: io_unit + + CALL check(ABS(actual - expected) < 1.0E-12_real_8, label, success, io_unit) + END SUBROUTINE check_real + +END MODULE dbcsr_test_symmetry diff --git a/tests/dbcsr_unittest1.F b/tests/dbcsr_unittest1.F index 18a39af19b1..776a5e2d6a6 100644 --- a/tests/dbcsr_unittest1.F +++ b/tests/dbcsr_unittest1.F @@ -25,6 +25,7 @@ PROGRAM dbcsr_unittest_1 mp_world_finalize, & mp_world_init, mp_comm_type USE dbcsr_test_add, ONLY: dbcsr_test_adds + USE dbcsr_test_symmetry, ONLY: dbcsr_test_hermitian_symmetry USE dbcsr_test_methods, ONLY: dbcsr_reset_randmat_seed USE dbcsr_test_multiply, ONLY: dbcsr_test_multiplies USE dbcsr_types, ONLY: dbcsr_mp_obj @@ -74,6 +75,8 @@ PROGRAM dbcsr_unittest_1 ! run tests + success = dbcsr_test_hermitian_symmetry(group, npdims, io_unit) .AND. success + ! add ----------------------------------------------------------------------- CALL dbcsr_test_adds("add_1", &