dbcsr_add_anytype Subroutine

   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, ndata, chunk, first, count
      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

         ! Equal, packed CSR layouts can be updated directly. This is entirely
         ! rank-local, so each rank may independently select the fast path.
         IF (same_data_layout(matrix_a, matrix_b)) THEN
            ndata = dbcsr_data_get_size_referenced(matrix_a%data_area)
            ! Reuse the existing type-specific AXPY operations without creating
            ! work matrices or rebuilding the CSR index.
!$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(chunk, first, count) &
!$OMP              SHARED(matrix_a, matrix_b, data_type_a, ndata, do_scale, my_beta_scalar) &
!$OMP              SCHEDULE(STATIC) IF(ndata > 131072)
            DO chunk = 0, (ndata - 1)/32768
               first = chunk*32768 + 1
               count = MIN(32768, ndata - first + 1)
               SELECT CASE (data_type_a)
               CASE (dbcsr_type_real_4)
                  CALL dbcsr_update_contiguous_blocks_s(matrix_a, matrix_b, first, first, count, &
                                                        do_scale, my_beta_scalar, .TRUE., 1)
               CASE (dbcsr_type_real_8)
                  CALL dbcsr_update_contiguous_blocks_d(matrix_a, matrix_b, first, first, count, &
                                                        do_scale, my_beta_scalar, .TRUE., 1)
               CASE (dbcsr_type_complex_4)
                  CALL dbcsr_update_contiguous_blocks_c(matrix_a, matrix_b, first, first, count, &
                                                        do_scale, my_beta_scalar, .TRUE., 1)
               CASE (dbcsr_type_complex_8)
                  CALL dbcsr_update_contiguous_blocks_z(matrix_a, matrix_b, first, first, count, &
                                                        do_scale, my_beta_scalar, .TRUE., 1)
               END SELECT
            END DO
!$OMP           END PARALLEL DO
            IF (PRESENT(flop)) flop = flop + 2_int_8*ndata
            CALL timestop(handle)
            RETURN
         END IF

         ! 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