outer_product_arrays Module Function

  module function outer_product_arrays(a, b) result(c)
    !! Outer product of two autodiff arrays
    implicit none
    class(array_type), intent(in), target :: a, b
    type(array_type), pointer :: c

    integer :: i, j, s

    ! check shapes
    if(size(a%shape).ne.1 .or. size(b%shape).ne.1)then
       call stop_program("dot_product_arrays: only 1D arrays supported")
    elseif(size(a%val,2).ne.size(b%val,2))then
       call stop_program("dot_product_arrays: array length mismatch")
    end if

    c => a%create_result(array_shape = [size(a%val,1), size(b%val,1), size(a%val,2)])
    ! outer product 1D array by using shape to swap dimensions
    do concurrent(s=1:size(a%val,2))
       do concurrent(i=1:size(a%val,1), j=1:size(b%val,1))
          c%val(i + (j-1)*size(a%val,1),s) = a%val(i,s) * b%val(j,s)
       end do
    end do

    c%get_partial_left => get_partial_outer_product_left
    c%get_partial_right => get_partial_outer_product_right
    c%is_sample_dependent = a%is_sample_dependent
    if(a%requires_grad .or. b%requires_grad) then
       c%requires_grad = .true.
       c%is_forward = a%is_forward .or. b%is_forward
       c%operation = 'outer_product'
       c%left_operand => a
       c%right_operand => b
       c%owns_left_operand = a%is_temporary
       c%owns_right_operand = b%is_temporary
    end if
  end function outer_product_arrays