dbcsr_dot_d Subroutine

      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