Solving sparse matrix on GPU and memory problems (even with free memory available)

Hi there,

I am running a backslash computation A\b where A is an array of type sparse and b is a vector, both stored as gpuArrays. It works fine for small matrices, but the following error is given for a matrix A of order 2e5:

_ Error using \ The GPU failed to allocate memory. To continue, reset the GPU by running 'gpuDevice(1)'. If this problem persists, partition your computations into smaller pieces._

KGR_gpu=gpuArray(KGR);
FR_gpu=gpuArray(FR);
sol=KGR_gpu\FR_gpu;

This same computation works fine when gpuArrays are not used. In fact, the size of A is only 50Mb and the GPU is a Nvidia GTX 1050, 4Gb. Following is the result of a gpuDevice call after this error:

   CUDADevice with properties:

                      Name: 'GeForce GTX 1050'
                     Index: 1
         ComputeCapability: '6.1'
            SupportsDouble: 1
             DriverVersion: 9.2000
            ToolkitVersion: 8
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [2.1475e+09 65535 65535]
                 SIMDWidth: 32
               TotalMemory: 4.2950e+09
           AvailableMemory: 3.3949e+09
       MultiprocessorCount: 5
              ClockRateKHz: 1493000
               ComputeMode: 'Default'
      GPUOverlapsTransfers: 1
    KernelExecutionTimeout: 1
          CanMapHostMemory: 1
           DeviceSupported: 1
            DeviceSelected: 1

It seems that this is an issue with the backslash operator since both A and b can be stored on the GPU without problems.

Any thoughts on how to solve large sparse systems on the GPU? Thanks!

Regards, Paulo

Answer 1

[Walter Roberson]

The result of \ between two sparse arrays is generally a dense array. For example sprand(1000,1000,.01) \ sprand(1000,1000,.01) gave me a result with fill fraction of 0.997031

[Paulo Ribeiro]

Dear Walter, thanks. But b is a vector. So it's a sparse array A and a vector b, both stored in the GPU. Perhaps this specific GPU is not able to perform these operations using double arrays? I know that single performance works great, but cant be used with sparse.

[Walter Roberson]

>> t = sprand(1000,1000,.01) \ sprand(1000,1,.01); nnz(t)./numel(t)
ans =
     1
>> whos t
  Name         Size            Bytes  Class     Attributes

    t         1000x1             16016  double    sparse

Notice this is twice the storage size that would be required for a non-sparse array with the same number of elements, due to the overhead of storing sparse arrays.

Answer 2

[Joss Knight]

Hi Paulo. I must admit, I'm not extremely familiar with the behaviour of the matrix factorization we use to implement the sparse direct solve; however it wouldn't surprise me if the result is quite dense. It might be interesting to (on the CPU) look at the density of the QR factors using the qr function, for your particular input. Certainly, when I did it on a random matrix with 10% fill, the Q factor was nearly completely dense and R factor was 50%.

>> A = sprand(1000, 1000, 0.1);
>> [Q,R] = qr(A);
>> nnz(Q)/numel(Q)
ans =
    0.9903
>> nnz(R)/numel(R)
ans =
    0.5005

For LU, both factors are 50% dense.

Obviously random sparse matrices don't properly reflect the structure of real sparse matrices, so your problem would be different. But it's not unreasonable to surmise that the intermediate factors might be very large.

To circumvent such problems a normal approach would be to use an iterative solver like gmres, bicg, pcg, cgs, lsqr etc. It is not uncommon for these to converge quicker than the direct solve can, especially if you can give them a good preconditioner.

[Paulo Ribeiro]

Thanks Joss,

Very nice. I am at this present time trying an iteractive solver combined with the GPU. Would you recommend any additional steps to this code?

L=ichol(KGR); 
sol= pcg(KGR_gpu,FR_gpu,1e-1,1e6,L,L');

For the time being this same code runs whenever matrix KGR and vector FR are not stored at the GPU. But for some reason, it does not work with gpuArrays:

Error using gpuArray/pcg (line 58) When the first input argument is a sparse matrix, the second preconditioner cannot be a matrix. Use functions for both preconditioners, or multiply the precondition matrices into one matrix.

By following the recommendation for a product of the preconditioners it follows:

_ _ _ Undefined function 'itermsg' for input arguments of type 'char'.

Error in gpuArray/pcg>iPcgBuiltinWrapper (line 119) itermsg('pcg',tol,maxit,0,flag,iter,NaN);

Error in gpuArray/pcg (line 63) [varargout{1:nargout}] = iPcgBuiltinWrapper(varargin{:});_____________

Thanks!

Regards,

Paulo

[Joss Knight]

That's a bug! Thanks for finding that. This is hit if the RHS vector is all zeros, which is a pathological case, so with a real RHS you'll avoid this.

[Paulo Ribeiro]

Hi Joss and many others that may follow this thread. The problem is solved (not yet - please check Joss comments) with preallocation on the GPU:

KGR_gpu=gpuArray(sparse(n,n))    % where n is the order of KGR_gpu
KGR_gpu=KGR;                     % KGR on the CPU is stored on GPU

With this workaround, the backslash and the iterative methods are working. Since I am using a Nvidia GTX 1050 (not suitable for double precision) the backslash with KGR_gpu and KGR are almost equivalent in time:

 tic
 KGR_gpu\FR
 toc

 tic
 KGR\FR
 toc

Where FR is a vector. Any experience with the same problem using double precision cards such as Nvidia Tesla? Are speedups expected with the backslash and KGR on the GPU?

Thanks!

[Joss Knight]

I'm glad you're making progress but it doesn't look as though you're doing anything on the GPU now. This code:

KGR_gpu=gpuArray(sparse(n,n))    % where n is the order of KGR_gpu
KGR_gpu=KGR;                     % KGR on the CPU is stored on GPU

may allocate a GPU array, but it then immediately frees it and overwrites the variable with the original KGR variable on the CPU. After this code, KGR_gpu and KGR are the same CPU array and the computation will happen on the CPU in both cases.

[Paulo Ribeiro]

Thanks Joss. Indeed this is a mistake. Made the same computations with:

KGR_gpu=gpuArray(sparse(length(KGR),length(KGR)));
KGR_gpu=gpuArray(KGR);

and the same error of the beginning of this thread arises once again:

Error using  \ The GPU failed to allocate memory. To continue, reset the GPU by running 'gpuDevice(1)'. If this problem
  persists, partition your computations into smaller pieces.

This is still very strange since KGR has only 0.20Gb and the GPU is a GTX 1050 with 4Gb. I will post updates of this same computation with an iterative solver.

Regards,

Paulo

[Paulo Ribeiro]

This same problem is solved without any issues with a pcg iterative solver:

sol=pcg(KGR_gpu,FR,1e-5,1e5);

and is faster than the CPU version:

sol=pcg(KGR,FR,1e-5,1e5);

For reference: GPU elapsed time = 47s; CPU elapsed time = 323s; Speedup = 6.9

[Joss Knight]

Good to know this is working out for you. If performance is an issue for you you should keep experimenting with the supported solvers: gmres, pcg, bicg, bicgstab, cgs, lsqr. Each will have different properties and one may converge faster for your problem than another.

Another trick is to pass your system matrix directly as the preconditioner:

pcg(KGR_gpu, FR, 1e-5, 1e5, KGR_gpu);

The GPU solvers (currently) use ILU to factorize the matrix and attempt to precondition the system. Sometimes this makes no difference but often it can dramatically reduce the number of iterations to convergence.

[Paulo Ribeiro]

Thanks Joss,

For this specific case the pcg performs better. No further improvement was achieved using a preconditioner. What still amazes me is the fact that a 0.20Gb matrix can't be solved using the backslash operator in a nominal 4Gb GPU. On the other hand, while monitoring the CPU version of this code (on Windows task manager), I can see that MATLAB bites almost 4Gb of RAM during the backslash operation. That might explain something.

I wonder how far can you get even with high-end GPU's.

Regards,

Paulo

[Walter Roberson]

I seem to recall that some memory is required to marshal the data on the GPU, which uses arrays in a different order than MATLAB uses. I do not recall the details at the moment, but potentially you might not be able to create an output array larger than half of your GPU memory.