Out of GPU memory

[albertomercurio]

Hello,

By running the following code I get the Out of GPU memory error.

using LinearAlgebra
using SparseArrays
using CUDA
using Krylov
using KrylovPreconditioners

N = 1500000
A = 10*CUSPARSE.CuSparseMatrixCSR(sprand(ComplexF64, N, N, 1/N)) + I
b = CUDA.CuArray(rand(ComplexF64, N))

Pℓ = Diagonal(CuVector(1 ./ norm.(eachcol(SparseMatrixCSC(A)))))
x, stats = gmres(A, b, N=Pℓ, ldiv=false)

I'm not surprised that this can happen, given the huge size of the matrix. But I'm surprised that this happens when the algorithm has already started and a few iterations have been already done. This also happens if I use the KrylovOperator, that should improve the performance.

Am I doing something wrong?

[amontoison]

@albertomercurio It is the main drawback of GMRES, it requires more storage at each iteration:
https://jso.dev/Krylov.jl/dev/storage/

So Krylov.jl allocates additional storage when the current number of iterations is greater than the memory parameter.

To avoid this issue, you need to use DQGMRES or try the restarted option of GMRES.

[albertomercurio]

Which one do you recommend?

[amontoison]

I suggest DQGMRES with as much "memory" as you can use (number of vectors in the Krylov basis).
You just to find the bigger k in https://jso.dev/Krylov.jl/dev/storage/ that allows to store the DqgmresWorkspace on your GPU.

[albertomercurio]

Thank you, I don't have the error anymore. I managed to set the memory up to 400.

However, the solver keeps running for hours, without finding a solution. Is there something we can do? Maybe other algorithms are better?

[amontoison]

For random linear systems, Krylov solvers don't perform well in practice.
What I recommend is (i) scaling, (ii) preconditioning, (iii) try bicgstab and restart + warm-start if you encounter a breakdown.

[albertomercurio]

Yes mine was a MWE. In my actual case I have a Liouvillian matrix, so not random. It is a non hermitian matrix, but it has one steady state that can be extracted by linear solving

[amontoison]

Is it possible to dump the matrix MatrixMarket format + right-hand side? (.mtx)
I can have a look.
Did you also try to use the ILU(0) preconditioner (from CUSOLVER) on your problem?
It is available as an operator in KrylovPreconditioners.jl.

[albertomercurio]

Sorry for my late reply.

Did you also try to use the ILU(0) preconditioner (from CUSOLVER) on your problem?

Yes I tried it as a right preconditioner.

I could send you the MatrixMarket format files, but I can send the simple code here. Tell me if you prefer to just have the matrix, but here you can play by increasing the matrix size by changing the Na and Nb parameters.

using LinearAlgebra
using SparseArrays
using QuantumToolbox
using CUDA
using LinearSolve

using QuantumToolbox: _dense_similar, _sparse_similar, allowed_setindex!

function generate_system(L)
    M_tmp = L.data
    N = prod(L.dimensions)
    weight = norm(M_tmp, 1) / N^2

    v0 = _dense_similar(M_tmp, N^2)
    fill!(v0, 0)
    allowed_setindex!(v0, weight, 1) # Because scalar indexing is not allowed on GPU arrays

    idx_range = collect(1:N)
    rows = _dense_similar(M_tmp, N)
    cols = _dense_similar(M_tmp, N)
    vals = _dense_similar(M_tmp, N)
    fill!(rows, 1)
    copyto!(cols, N .* (idx_range .- 1) .+ idx_range)
    fill!(vals, weight)
    Tn = _sparse_similar(M_tmp, rows, cols, vals, N^2, N^2)
    M_tmp = M_tmp + Tn

    return M_tmp, v0
end

# %%

L = 1f0

κ  = 0.14f0
γ  = 0.001f0 # 2.5e-8
ωb = 1f0
g0  = 0.36f0 / √L
Δ  = -1f0
F = 0.1f0 * √L

tlist = range(0, 4.0 / γ, 2000)

# %%

Na = 20
Nb = 20

a = tensor(destroy(Na), qeye(Nb)) |> CUSPARSE.CuSparseMatrixCSR{ComplexF32}
b = tensor(qeye(Na), destroy(Nb)) |> CUSPARSE.CuSparseMatrixCSR{ComplexF32}

# %%
# --------- FULL QUANTUM ----------

H = -Δ*a'*a + ωb*b'*b + g0*a'*a*(b + b') + F*(a + a')
c_ops = [√κ*a, √γ*b]
L0 = liouvillian(H, c_ops)

generate_system(L0)

The Liouvillian is entirely described by the L0 object, you can access the matrix with L0.data. The steadystate is found by solving

$$ \frac{d}{dt} \rho = L \rho = 0 $$

but numerically you get the trivial solution $\rho = 0$. That's why I define the generate_system, that solves

$$ (L + T) \rho = \mu \qquad \text{where} \qquad \mu = \begin{pmatrix}w \ 0 \ \dots \ 0\end{pmatrix}^T $$

with $w$ some weight.

Sorry for being too specific, just wanted to explain the context.