Quasi-Newton operators compatible with GPUs

I am wondering how much work it would be to make QN operators compatible with GPUs (CuArray for instance)?

Typically, my use case would be something like this
```
using CUDA, NLPModels, NLPModelsModifiers, NLPModelsTest
V = CuArray{Float64}
nlp = NLSLC(V)

CUDA.allowscalar()
list_QN = [LBFGSModel, LSR1Model, DiagonalPSBModel, DiagonalAndreiModel, SpectralGradientModel]
lnlp = list[1](nlp)
x = nlp.meta.x0
v = copy(x)
Hv = similar(x)
hprod!(lnlp, x, v, Hv)
```
These model modifiers internally call the QN operators from LinearOperators.jl, e.g. `op = LBFGSOperator(T, nlp.meta.nvar)`.