Dot product of DBCSR matrices
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(dbcsr_type), | intent(in) | :: | matrix_a |
DBCSR matrices DBCSR matrices |
||
type(dbcsr_type), | intent(in) | :: | matrix_b |
DBCSR matrices DBCSR matrices |
||
complex(kind=real_4), | intent(inout) | :: | trace |
the trace of the product of the matrices |
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