Hi all,

I am trying to apply an array of ITensors to an MPO exactly (without truncation). Here is a MWE.

```
using LinearAlgebra
using ITensors
num_sites=10
sites = siteinds("S=1/2", num_sites)
psi=randomMPS(sites,linkdims=23)
A=randomITensor(sites[1]',sites[2]',sites[1],sites[2])
B=randomITensor(sites[3]',sites[4]',sites[3],sites[4])
itarray=[A,B]
psi=apply(itarray,psi,method="naive",truncate=false)
```

This throws the following error.

```
MethodError: no method matching factorize(::ITensor, ::Vector{Index{Int64}}; method="naive", truncate=false, tags="Link,n=1", ortho="left")
Closest candidates are:
factorize(::ITensor, ::Any...; mindim, maxdim, cutoff, ortho, tags, plev, which_decomp, eigen_perturbation, svd_alg, use_absolute_cutoff, use_relative_cutoff, min_blockdim, singular_values!, dir) at ~/.julia/packages/ITensors/MnaxI/src/tensor_operations/matrix_decomposition.jl:727 got unsupported keyword arguments "method", "truncate"
factorize(::ITensor; kwargs...) at ~/.julia/packages/ITensors/MnaxI/src/tensor_operations/matrix_decomposition.jl:682
```

I have checked the documentation that this works for application of MPO to MPS. But is there an implementation of the same for an array of ITensors? It seems to work for an array of ITensors in v0.3.39 but doesn’t work on v0.3.57.

Also a follow-up question. Is there an in-place version of apply (e.g. apply!). Does `psi=apply(itarray,psi..)`

cause unnecessary allocation that can be avoided?

Thank you so much for this awesome library!