matmul_real2d Module Function

  module function matmul_real2d(a, b) result(c)
    !! Matrix multiplication of a real array and an autodiff array
    !! Computes C = a * b where a is autodiff (vector per sample) and b is a
    !! real 2D matrix. Equivalent to C(:,s) = b^T * a(:,s) for each sample s.
    implicit none
    class(array_type), intent(in), target :: a
    real(real32), dimension(:,:), intent(in) :: b
    type(array_type), pointer :: c
    type(array_type), pointer :: b_array

    integer :: s, i, rows, cols, num_samples

    rows = size(b, 1)
    cols = size(b, 2)
    num_samples = size(a%val, 2)

    c => a%create_result(array_shape = [cols, num_samples])
#ifdef USE_BLAS
    ! C(cols, S) = b^T(cols, rows) * a%val(rows, S)
    ! sgemm: transa='T' transposes b in-place, no temporary needed
    ! m = cols, n = S, k = rows
    ! lda = rows (leading dim of b before transpose), ldb = rows, ldc = cols
    call sgemm('T', 'N', cols, num_samples, rows, &
         1.0_real32, b, rows, a%val, rows, 0.0_real32, c%val, cols)
#else
    c%val = matmul(transpose(b), a%val)
#endif

    c%is_sample_dependent = a%is_sample_dependent
    c%get_partial_left => get_partial_matmul_left
    c%get_partial_right => get_partial_matmul_right
    c%get_partial_left_val => get_partial_matmul_left_val
    c%get_partial_right_val => get_partial_matmul_right_val
    c%get_partial_left_val_sum => get_partial_matmul_left_val_sum
    c%get_partial_right_val_sum => get_partial_matmul_right_val_sum
    if(a%requires_grad) then
       c%requires_grad = .true.
       c%is_forward = a%is_forward
       c%operation = 'matmul_scalar'
       c%left_operand => a
       c%owns_left_operand = a%is_temporary
    end if
    allocate(b_array)
    b_array%is_sample_dependent = .false.
    b_array%shape = shape(b)
    b_array%requires_grad = .false.
    call b_array%allocate(array_shape=[size(b,1), size(b,2), 1])
    b_array%val(:,1) = reshape(b, [size(b,1)*size(b,2)])
    c%right_operand => b_array
    c%owns_right_operand = .true.
  end function matmul_real2d