contract_tensors Subroutine

pure subroutine contract_tensors(A, B, C)

    real(dp), dimension(:), intent(in) :: A
    real(dp), dimension(:), intent(in) :: B
    real(dp), dimension(:), intent(inout) :: C

    integer :: x_a, y_a, z_a
    integer :: x_b, y_b, z_b
    integer :: x_c, y_c, z_c
    integer :: idx_a, idx_b, idx_c
    integer :: rank_a, rank_b, rank_c
    integer :: error_code

    rank_a = tensor_rank(A)
    rank_b = tensor_rank(B)
    rank_c = tensor_rank(C)

    error_code = 0

    if (rank_a < rank_b) then
        ! TODO: return control to caller with error code
        ! 'MiMiC ERROR: rank(A) must be equal to or larger than rank(B)'
        error_code = error_code + 1
    end if

    if (rank_c /= rank_a - rank_b) then
        ! TODO: return control to caller with error code
        ! 'MiMiC ERROR: rank(C) must be equal to rank(A) - rank(B)'
        error_code = error_code + 1
    end if

    do x_c = rank_c, 0, -1
        do y_c = rank_c - x_c, 0, -1
            z_c = rank_c - x_c - y_c
            idx_c = convert_index(x_c, y_c, z_c)
            do x_a = rank_a, x_c, -1
                do y_a = rank_a - x_a, y_c, -1
                    z_a = rank_a - x_a - y_a
                    if (z_a < z_c) cycle
                    idx_a = convert_index(x_a, y_a, z_a)
                    x_b = x_a - x_c
                    y_b = y_a - y_c
                    z_b = z_a - z_c
                    idx_b = convert_index(x_b, y_b, z_b)
                    C(idx_c) = C(idx_c) + trinomial(x_b, y_b, z_b) * A(idx_a) * B(idx_b)
                end do
            end do
        end do
    end do

end subroutine contract_tensors