dbcsr_hadamard_product Subroutine

   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