Matrix and vector multiplication of size using a CPU is very slow. Using GPU is much quicker but I need a way around the size limitation.

This is an extract from the core part of a project I am working on. It is essentially a Finite Difference implementation of Crank Nicholson method. If I run a non-GPU version app.U0 can easily be a vecor of size 64k. The matrices are of similar dimension. The non-gpu version can take a day to run on my overclocked AMD Ryzen 7 5700X which has 128GB memory. I thought I would try and use the GPU

>> gpuDevice
ans = 
  CUDADevice with properties:
                      Name: 'NVIDIA GeForce RTX 3060'
                     Index: 1
         ComputeCapability: '8.6'
            SupportsDouble: 1
     GraphicsDriverVersion: '536.99'
               DriverModel: 'WDDM'
            ToolkitVersion: 11.8000
        MaxThreadsPerBlock: 1024
          MaxShmemPerBlock: 49152 (49.15 KB)
        MaxThreadBlockSize: [1024 1024 64]
               MaxGridSize: [2.1475e+09 65535 65535]
                 SIMDWidth: 32
               TotalMemory: 12884377600 (12.88 GB)

I have little idea of the significance of this other than the 12+GB which I thought suggested I had lots of room. The test I ran using both on and off GPU showed a perfect match. The only difference being the CPU version for the small app.U0 took 50 seconds and the GPU version took 2!! Total win... except... when the xMesh reaches >53 the array contians only NaN except for the initial input data.

The LHS diagram shows why I need larger than 100x2x50 start vector.The rays emanating from the start and ricocheing off the boundary are a real effect and one that needs to be avoided. This means larger distance between boundaries. Apparently fine if I use the CPU but not on the GPU. Using stochastic data for the noise I might want 30+ runs and that means 30+days just to get the data.

                A = diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);
                sparse(A);
                ul = gpuArray(complex(ul));
                ur = gpuArray(complex((ur)));
                app.u0 = addNoise(app,startAgain,loop);%first line of array eventually
                U = gpuArray(complex(zeros(M+2,app.NJ+1))); %U = zeros(M+2,app.NJ); %output file
                U0 = gpuArray(complex(app.u0(2:end-1,:)));
                U1 = gpuArray(complex(size(U0)));
                UC = gpuArray(complex(size(U0)));
                D = app.dt/(2*app.h^2); % D 
                Bdy(1,:)= D*ul;
                Bdy(end,:)=D*ur;
                %Bdy = sparse(Bdy);
                Bdy = gpuArray(complex(Bdy));
                AB = gpuArray(complex(D*A));
                AA = sparse(1i*eye(size(A))+AB); AA = gpuArray(complex(AA));
                sdg = gpuArray(s*app.dt*app.gamma);%this is just a constant
                    for j=1:app.NJ
                        % This is the first shot at Crank Nicholson
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0...
                            - Bdy(:,j)-Bdy(:,j+1));
                        %which becomes a predictor corrector by a second
                        %run through
                        UC = (U1+U0)/2;
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC...
                            - Bdy(:,j)-Bdy(:,j+1));
                        U0=U1; U(:,j+1) = [ul(j+1);U1;ur(j+1)];
                    end

                disp(f); disp(showDateAndTime(app,2));           %quick indicator that a loop has finished
                %put start on data
                %U(:,1) = app.u0;%U = [app.u0 U];
                U = gather(U);
                U(:,1) = app.u0;

The 'for' loop is the heart of the matter. Is there any way around this issue without losing all the benefits of the GPU speed?

[Matt J]

It's not clear (to me) what you mean by the "issue" or the "size limitation". Is there some difference you see in the GPU and CPU versions other than execution speed? Are you getting out of memory error messages?

[Jonathan Wharrier]

The answers I get running on the CPU and the GPU are absolutely identical until the size of the start has a dimension of >53 which means that U0 start vector is 2x53x100 by 1 in dimension at which point I get returned NaN in all locations. The CPU has no such problem and I have run up to xMesh +/-320 (64k by 1). The x100 is because the space slice for the Crank Nicholson is 0.01. There seems to be a size issue that I do not understand. Less than ~10600 for the dimension of the input on the GPU and it is fine.

When MATLAB has finally reloaded on this machine (Tx MS!?) I can find out if it does this to the time dimension. It is also discretised. That may make things clearer.

[Bruno Luong]

Can you provide a MWE where GPU returns NaN?

[Jonathan Wharrier]

I have attached an output file which is U in the code. This one is 10801x501 complex double. Column 1 is U0. Rows 1 and 10801 are zero (ul & ur) because of the requirement for compact support. The rest is where the GPU has put U1 for each increment of time. The result is output as a .mat file and the 3-D picture (#2 above) is the outcome. There is no error recorded (!).

    load('U_L54run1P_1.mat');

The only difference between runs is that XRange for the sample was 54. Using an xRange of 53 gave the complete 3D above.

[Bruno Luong]

I need the INPUTS AA , AB, U0, U1, sdg, Bdy in order to run the for-loop bellow. I want to see when and why it returns NaN. You can give them in standard array so I can run on CPU. I can cast my selft to run on GPU.

for j=1:app.NJ
    % This is the first shot at Crank Nicholson
    U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0...
        - Bdy(:,j)-Bdy(:,j+1));
    %which becomes a predictor corrector by a second
    %run through
    UC = (U1+U0)/2;
    U1=AA\(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC...
        - Bdy(:,j)-Bdy(:,j+1));
    U0=U1; U(:,j+1) = [ul(j+1);U1;ur(j+1)];
end

The output colapsing to NaN is not terribly useful.

[Joss Knight]

Do not use sparse backslash on a gpuArray, it is nearly always slower than CPU. Either use CPU or use a iterative solver like GMRES.

[Bruno Luong]

@Jonathan Wharrier here is why nobody can effectively help you:

                A = diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);

Unrecognized function or variable 'M'.

                sparse(A);
                ul = gpuArray(complex(ul));
                ur = gpuArray(complex((ur)));
                app.u0 = addNoise(app,startAgain,loop);%first line of array eventually
                U = gpuArray(complex(zeros(M+2,app.NJ+1))); %U = zeros(M+2,app.NJ); %output file
                U0 = gpuArray(complex(app.u0(2:end-1,:)));
                U1 = gpuArray(complex(size(U0)));
                UC = gpuArray(complex(size(U0)));
                D = app.dt/(2*app.h^2); % D 
                Bdy(1,:)= D*ul;
                Bdy(end,:)=D*ur;
                %Bdy = sparse(Bdy);
                Bdy = gpuArray(complex(Bdy));
                AB = gpuArray(complex(D*A));
                AA = sparse(1i*eye(size(A))+AB); AA = gpuArray(complex(AA));
                sdg = gpuArray(s*app.dt*app.gamma);%this is just a constant
                    for j=1:app.NJ
                        % This is the first shot at Crank Nicholson
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0...
                            - Bdy(:,j)-Bdy(:,j+1));
                        %which becomes a predictor corrector by a second
                        %run through
                        UC = (U1+U0)/2;
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC...
                            - Bdy(:,j)-Bdy(:,j+1));
                        U0=U1; U(:,j+1) = [ul(j+1);U1;ur(j+1)];
                    end

                disp(f); disp(showDateAndTime(app,2));           %quick indicator that a loop has finished
                %put start on data
                %U(:,1) = app.u0;%U = [app.u0 U];
                U = gather(U);
                U(:,1) = app.u0;

[Jonathan Wharrier]

I have cut the code out of the app and added enough data so that the loop runs. Changing xRange from 50 to 54 reproduces the problem. The output <U> becomes mainly NaN the non NaN are the boundary values (all zero) and the start value.

           % CRANK NICHOLSON VERSION
            xRange = 50; %This is the critical value which defines the size of the mesh grid
            % crashes at 54
            dx = 0.01; %mesh grid increment size
            x = -xRange:dx:xRange;%mesh grid
            T = 5; %time for run
            dt = 0.01; %Time increment size
            NJ=T/dt; %number of iterations
            t= dt*(0:NJ); %time vector
            M=length(x);    % length of the mesh
            M=M-2;          % active length minus end points
                                
            S  = @(ex) 1*sech(sqrt(1/2)*ex).*exp(1i*2*ex); %creates Soliton
            u0 = S(x)'; %Start data - note transform
            s=1; %constant
            ul = zeros(size(t)); %set boundary values
            ur = zeros(size(t)); %boundaries are zero
                reset(gpuDevice);
                Bdy=zeros(M,NJ+1);
                A = diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);
                sparse(A);
                U = gpuArray(complex(zeros(M+2,NJ+1))); %U = zeros(M+2,app.NJ); %output file
                U0 = gpuArray(complex(u0(2:end-1,:)));
                UC = gpuArray(complex(size(U0)));
                D = dt/(2*dx^2); % D 
                Bdy(1,:)= D*ul;
                Bdy(end,:)=D*ur;
                %Bdy = sparse(Bdy);
                url = gpuArray(complex(ul));
                urr = gpuArray(complex((ur)));
                Bdy = gpuArray(complex(Bdy));
                AB = gpuArray(complex(D*A));
                AA = sparse(1i*eye(size(A))+AB); AA = gpuArray(complex(AA));
                sdg = gpuArray(1*dt*1);%this is just a constant
                    for j=1:NJ
                        % This is the first shot at Crank Nicholson
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0...
                            - Bdy(:,j)-Bdy(:,j+1));
                        %which becomes a predictor corrector by a second
                        %run through
                        UC = (U1+U0)/2;
                        U1=AA\(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC...
                            - Bdy(:,j)-Bdy(:,j+1));
                        U0=U1;
                        %U(:,j+1) = [ul(j+1);U1;ur(j+1)];
                        U(:,j+1) = [url(j+1);gather(U1);urr(j+1)];
                        
                    end
                U = gather(U);
                U(:,1) = u0;
            close all;
            surf(abs(U));

[Jonathan Wharrier]

The run which took 52 minutes on the CPU took less than 2 on the GPU even using the backslash solver.

[Joss Knight]

Yes. It looks like AA is tridiagonal which is a special case and so runs fast.

Answer 1

[Bruno Luong]

I run your (slighly modified) code with xRange = 100 (EDIT) and get the finite result

My config: R2023a, Windows 11, Laptop 32 Gbytes, Intel i9 12th generation and Nvidia RTX 3060.

I would look for updateing your MATLAB, graphic card drivers, or eventually HW incompatibility issue.

% CRANK NICHOLSON VERSION
xRange = 100; %This is the critical value which defines the size of the mesh grid
% crashes at 54
dx = 0.01; %mesh grid increment size
x = -xRange:dx:xRange;%mesh grid
T = 5; %time for run
dt = 0.01; %Time increment size
NJ=T/dt; %number of iterations
t= dt*(0:NJ); %time vector
M=length(x);    % length of the mesh
M=M-2;          % active length minus end points

S  = @(ex) 1*sech(sqrt(1/2)*ex).*exp(1i*2*ex); %creates Soliton
u0 = S(x)'; %Start data - note transform
s=1; %constant
ul = zeros(size(t)); %set boundary values
ur = zeros(size(t)); %boundaries are zero
reset(gpuDevice);
Bdy=zeros(M,NJ+1);
A = spdiags(ones(3*M,1)*[1 -2 1],[-1 0 1],M,M); % modified by Bruno %diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);
U = gpuArray(complex(zeros(M+2,NJ+1))); %U = zeros(M+2,app.NJ); %output file
U0 = gpuArray(complex(u0(2:end-1,:)));
UC = gpuArray(complex(size(U0)));
D = dt/(2*dx^2); % D
Bdy(1,:)= D*ul;
Bdy(end,:)=D*ur;
%Bdy = sparse(Bdy);
url = gpuArray(complex(ul));
urr = gpuArray(complex((ur)));
Bdy = gpuArray(complex(Bdy));
AB = gpuArray(complex(D*A));
AA = 1i*speye(size(A))+AB; 
AA = gpuArray(AA); % modified by Bruno
sdg = gpuArray(1*dt*1);%this is just a constant
for j=1:NJ
    % This is the first shot at Crank Nicholson
    U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0...
        - Bdy(:,j)-Bdy(:,j+1));
    %which becomes a predictor corrector by a second
    %run through
    UC = (U1+U0)/2;
    U1=AA\(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC...
        - Bdy(:,j)-Bdy(:,j+1));
    U0=U1;
    %U(:,j+1) = [ul(j+1);U1;ur(j+1)];
    U(:,j+1) = [url(j+1);gather(U1);urr(j+1)];

end
U = gather(U);
U(:,1) = u0;
close all;
surf(abs(U),'EdgeColor','none');

[Bruno Luong]

BTW the same code (attached) takes

Elapsed time is 2.059563 seconds.

on my CPU.

So I don't get anything like you describe.

[Jonathan Wharrier]

I am to say the least somewhat puzzled but nonetheless grateful for your efforts.

My system is

AMD Ryzen 7 5700X 8-Core Processor 3.40 GHz

128GB RAM

Graphics: 'NVIDIA GeForce RTX 3060'

If I run 50 in the app it takes 2s but as the stand alone it is taking nearer 8s.

When I put in xRange=100 I get a different result from xRange = 54. the system throws an error. I am watching Task manager CPU as they run. The 54 runs do not seem to have a problem with the GPU whereas at 100 one can understand why the out of memory as I watch the peak value hit the top. So why does yours not fail at this point I wonder.

Error using  + 
Out of memory on device. To view more detail about available memory on the GPU, use 'gpuDevice()'. If the problem
persists, reset the GPU by calling 'gpuDevice(1)'.

Error in GPUTest (line 38)
AA = sparse(1i*eye(size(A))+AB); AA = gpuArray(complex(AA));

[Jonathan Wharrier]

I think that you arer right in that I may have a hardware problem. I have just seen your output and realised that mine is faulty. This is my surf. figure

[Jonathan Wharrier]

One final question. You have an RTX 3060. What is the driver you are using ATM? I'd like to check mine against a working one before I get on to the company that supplied mine.

[Bruno Luong]

GPU time for xRange = 100

>> CPUtest
Elapsed time is 0.871546 seconds.

[Bruno Luong]

I recall there are some discussion threads about MATLAB performance with AMD processors.

As I have Intel I never pay great attention to them.

[Jonathan Wharrier]

FYI according to Google AI the problem with MATLAB and AMD were apparently solved for 2020a. This is clearly an NVIDIA problem so I'll bother them! Thanks for the support.

[Bruno Luong]

But your Nvidia card and mine looks like the same (Mine is Gforce RTX 3060 Laptop GPU). Do you have the latest driver?

Not sure what else that could affect the gpu calculation.

https://www.mathworks.com/help/parallel-computing/gpu-computing-requirements.html

[Jonathan Wharrier]

Mines not the laptop but it is an RTX Geoforce 3060...

I am really puzzled. I have looked at each of the Gpu arrays and everything is fine until the first attempt to solve the equation with the line

    U1=AA\(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0);%...
    %- Bdy(:,j)-Bdy(:,j+1));

I have replaced the boundary with zeros and simply deleted Bdy for the same reason. Same thing. All the values in AA; U0;AB have the right number of non zero componets but as 53->54 we get a fail in the solution and U1 reverts to NaN but yours does not!! I have looked at NVIDIA and loaded the very latest drivers, both the gaming one and the arty (!) one and both do the same and yet yours works. that leaves some issue with AMD possibly? I am at a loss. Clearly I can use the CPU but that is back to a sigle run taking 20+ hours and for a stochastic examination I need 30+ runs!? Not feasible.

[Bruno Luong]

I don't know if it's matter but I prefer not using those GameReady driver, rather using the stable one approved by MS.

My driver is 31.0.15.2904 dated from 19 may 2023

[Jonathan Wharrier]

I did a rebuid over the weekend, nothing to do with this but MS loaded a driver and I had the same problem. It is why I tried the other two drivers. It has to connect with the use of the \ command to solve the equation. Everything on the RHS appears to be correct and yet with xRange increased by one the answer ends up NaN. I am wondering if it is a strange issue with the graphics card because if it were a problem with "\" then why would you not have it?...

[Bruno Luong]

You should instrumental the for-loop check when this happens

for ...
    sol = AA \ rhs;
    if any(isnan(sol))
        save debug.mat AA rhs
        return
    end
end

Note there are two "\" in the loop.

Then gives that file debug.mat to TMW suport for them to help you for deeper diagnostic.

[Jonathan Wharrier]

Who are TMW? I am not familiar with this abbreviation. However I ran your code suggestion. The result is fascinating even if I have no idea why it happens.I inserted

if any(isnan(U1))
    save debug.mat AA AB U0;
    return
end

and found that it fell over with xRange = 54 as before. However I now have a copy of a debug file (attached)with the vectors/matrices in it. Examinig it is quite difficult but I found that surf(abs(U),'EdgeColor','none'); is very revealing. This is the result for AA, AB and U0.

Which are as one would expect. The first two are diagonal matricesand the third is the Solon starter.It failed on the very first attempt to calculate the result. The result was simply a vector of NaN. So the RHS is OK but the calculation is not.

[Bruno Luong]

You debug.mat is not useful for debugging, what we (the person who makes disgnostic) need is rhs in debug.mat

rhs = (1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC - Bdy(:,j)-Bdy(:,j+1));
sol = AA\rhs; 
    if any(isnan(sol))
        save debug.mat AA rhs
        return
    end

Save the INPUTS of the operation that creates an NaN output. If you save the NaN output we can only observe the damage.

And I guess it creats NaN on your computer and not mine. So you have to submlit a bug report with that file and all the details about your computer HW.

[Jonathan Wharrier]

OK. My apologies for misunderstanding. Put it down to my venerable age :)

rhs is exactly as expected and so is AA. Here is the interesting thing. The solution calculated as U1 is all NaN BUT!! if I load the debug file and and type in the command window

q = AA\rhs; plot(abs(q));

I get a sensible answer that plots to a Soliton...

rhs = (1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0);
U1=AA\rhs;
%(1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0);%...
%- Bdy(:,j)-Bdy(:,j+1));
if any(isnan(U1))
    save debug.mat U1 AA AB U0 rhs;
    return
end

[Bruno Luong]

It sounds like a buffer overrun, a memory issue or collision. You really need someone who can investigate.

BTW why you put gather in the last statement of the for-loop

U(:,j+1) = [url(j+1);gather(U1);urr(j+1)];

Can you leave U1 on the gpu side?

[Jonathan Wharrier]

because U is not on the gpu however it is unecessary ...

I even tried keeping U1 on the CPU side but it made no difference

So I am left with rhs is correctly calculated and AA is also correct bu once I pass a certain size (xRange > 53) the expression AA\rhs produces NaN on my system but not on yours. What I have done is the following...

Firstly I did a new run picking arbitrarily xRange =60. It failed. Double clicking on U1 produces a list of NaN. I cleared the workspace then ran the following by habnd in the command window...

1. load('debug.mat') this loads the debug file and AA; AB; rhs; U0 & U1
2. surf(abs(AA),'EdgeColor','none') provides graph of tri diagonal matrix as before.
3. surf(abs(AB),'EdgeColor','none') provides graph of matrix as before.
4. plot(abs(U0)) provides graph of matrix as before.
5. plot(abs(rhs)) provides graph of matrix as before.
6. issparse(AA) logical 1 (yes), it does not say this in workspace, just gpuArray.
7. issparse(AB) logical 1 (yes), it does not say this in workspace, just gpuArray.
8. sol = AA\rhs; this failed in the program for U1.
9. plot(abs(sol)) produces non NaN and shows a graph of a slightly moved soliton as required.
10. isgpuarray(sol) logical 1 (yes),

Apparently I can get the gpu to do this by hand one at a time! I am at a total loss. The following code runs if xRange <= 53 and provides a graph of a Soliton translated by 20 units and a surf diagram. Values of 54 and above fail for no discernible reason that I can fathom though the calc above (by hand as it were) seems to provide the correct answer.

           tic
            % CRANK NICHOLSON VERSION
            xRange = 50; %This is the critical value which defines the size of the mesh grid
            % crashes at 54
            dx = 0.01; %mesh grid increment size
            x = -xRange:dx:xRange;%mesh grid
            T = 5; %time for run
            dt = 0.01; %Time increment size
            NJ=T/dt; %number of iterations
            t= dt*(0:NJ); %time vector
            M=length(x);    % length of the mesh
            M=M-2;          % active length minus end points
                                
            S  = @(ex) 1*sech(sqrt(1/2)*ex).*exp(-1i*2*ex); %creates Soliton
            u0 = S(x)'; %Start data - note transform
            s=1; %constant
                reset(gpuDevice);
                A = diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);
                A = sparse(A);
                U = gpuArray(complex(zeros(M+2,NJ+1)));
                U0 = gpuArray(complex(u0(2:end-1,:)));
                UC = gpuArray(complex(size(U0)));
                D = dt/(2*dx^2); % D 
                DA = sparse(D*A);
                AB = sparse(gpuArray(complex(DA)));
                AA = sparse((1i*eye(size(A))+AB));
                AA = gpuArray(complex(AA));
                
                sdg = gpuArray(1*dt*1);%this is just a constant
                    for j=1:NJ
                        % This is the first shot at Crank Nicholson
                        rhs = (1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0);
                        U1=AA\rhs;
                            if any(isnan(U1))
                                save debug.mat U1 AA AB U0 rhs;
                                return
                            end
                  %which becomes a predictor corrector by a second
                  %run through
                        UC = (U1+U0)/2;
                        rhs=(1i.*U0 - AB*U0 - sdg*(UC.*conj(UC)).*UC); %...
                        U1=(AA\rhs);
                        U0=U1;
                        U(:,j+1) = [0;U1;0];
                    end
               U(:,1) = u0;
            toc
            figure(1);
            surf(abs(U),'EdgeColor','none');
            figure(2)
            plot(x,abs(u0));
            grid on; hold on;

[Joss Knight]

This could well be a bug in NVIDIA sparse direct solver that MATLAB uses, specific to your hardware (either the 3060 or the specific limited environment of a laptop with power and memory constraints). We can try to reproduce it on another 30xx card but I suspect this will also require the WDDM (Windows) driver.

What is it about xRange that is truly important for you? Almost certainly, if there is an environment-related issue it will be due to the size of your matrix, not the values, so can you just use a different dx grid increment for the larger values so your matrices stay the same size?

I have also confirmed that you get the same answers for iterative solvers if you use AA as a preconditioner. This is because the GPU iterative solvers use ILU for preconditioning and ILU gives an exact answer for tridiagonal matrices. So if you substitute, for instance, lsqr(AA,rhs,[],[],AA) for AA\rhs you should get the same answer. I also found that gpuArray(gather(AA)\gather(rhs)) was faster still. I agree that none are as fast as the original direct solve on the GPU, but I am just trying to find a way to unblock you. I used the CPU option on your script with xRange=54 and it went from 2.6s to 4s, so slower, but still pretty fast.

[Jonathan Wharrier]

I have got it working after a fashion. By removing ALL sparse matrices and setting T = 0.01 (one iteration) it now works for 'all' values of xRange. However it is completely useless.

>> GPUTest (xR = 53)

Elapsed time is 37.947680 seconds.

>> GPUTest (xR = 60)

Elapsed time is 55.617800 seconds.

>> GPUTest (xR = 100)

Elapsed time is 246.843792 seconds.

[Jonathan Wharrier]

My computer is a desktop

AMD Ryzen 7 5700X 8-Core Processor 3.40GHz 128GB RAM

CUDADevice with properties:

Name: 'NVIDIA GeForce RTX 3060'
Index: 1
ComputeCapability: '8.6'
SupportsDouble: 1
GraphicsDriverVersion: '536.99'
DriverModel: 'WDDM'
ToolkitVersion: 11.8000

xRange is crucial. If you look at the first 3D graph you can see reflections from the sides. I have run using xRange = 320 !! It took a day on a CPU for one. I need, for stochastic purposes, 30+ runs. You can see why this is a problem. Basically I need a scheme to run through Crank-Nicholson twice as a predictor corrector and then do this for the equivalent of about 25 seconds (2500 iteration). Also 0.01 for dx is about as large as I can go as the errors start to climb.

I have noted your comment about solving the system with lsqr so I will give that a go.

TX

[Joss Knight]

Like I say, you are better off using the CPU direct solve for just that one calculation, that will give the fastest result, so gpuArray(gather(AA)\gather(rhs)). Certainly you will get slower results if you move everything to dense matrices,

I have run your code on an RTX3070 without any problems, but that was on Linux. I suspect I will need both the Ampere card and the WDDM driver to reproduce, but I don't have access to both.

[Bruno Luong]

@Joss Knight "This could well be a bug in NVIDIA sparse direct solver that MATLAB uses, specific to your hardware (either the 3060 or the specific limited environment of a laptop with power and memory constraints). We can try to reproduce it on another 30xx card but I suspect this will also require the WDDM (Windows) driver."

Note that I have the Nvidia RTX 3060 card on Laptop Windows 11 but the same gpu code works on my side.

[Jonathan Wharrier]

I'm old and my brain hurts so I'm going for a kip! I will use the advice provided to see what I can do in terms of size and speed. and let you know tomorrow.

Tx to both of you for the support.

[Jonathan Wharrier]

This morning has been very productive. Firstly everything ran fine up to 200 xRange. It took about 3 minutes. I tried to be clever and went for 240 and it fell over but I have spent a bit of time thinking about sparse. It transpires that the whilst the graphics card works sparse OK it appears to me that it manages memory poorly. The "out of memory" seemed to happen because the graphics card maintained both versions. I have been watching carefully the Task manager as the computer handles the code. I have therefore used the main memory to create the large matrices and to sparse them and then loaded them to the card. Turns out this works a treat. I now can get to xRange = 360. When I tried for 400 because the card appeared to be less than "full" it fell over because my CPU/RAm at 128GB ran out of memory!, not the card!. The time when I used the CPU only before I had the NVIDIA (xRange = 320) was 23+ hours for a single run whereas this morning 360 xRangetook 5 minutes. I am still puzzled as to why

U1 = gather(AA)\gather(rhs);

works though !?

At ~30s on the CPU is creatng matrices and I would presume that MATLAB has set up storage on the Graphics card. The CPU is doing the work.

2 minutes in the CPU has loaded the graphics card (GPU spike) and CPU memory usage drops from 128GB.

The GPU is not over extended in terms of memory use and runs for about 2 minutes

the last minutes are the CPU saving some very large files to the HD

[Bruno Luong]

" U1 = gather(AA)\gather(rhs);

works though !?

At ~30s on the CPU is creatng matrices and I would presume that MATLAB has set up storage on the Graphics card. The CPU is doing the work."

I'm lost, what is the point of using gpuArray at all then?

At one point you said cpu tooks much longer than gpu. I find it also strange because it did not find similar figure omy side.

What I'm retain from the discussion s this : Using sparse direct method (backslash then) gpuArray tooks half of the time of cpu, but not more (my experiens). However gpuArray fails (returns NaN after some size anditeration) on your side not on mine.

[Jonathan Wharrier]

My limited understanding of implementation suggests that the problem is AS. This is where, if xRange (300+) is very large that the program runs out of memory. However, trying to create A or AD or AS on the graphics card runs out of memory much faster (54+). However, this all sounds fine until I realised that U, which is the output, sits quite happily in the gpu card memory. It is the only complete matrix. All the others are well zeroed and therefore sparse or else they are simply vectors. So everything is now on the graphics card after a lot of one off calculations on setup. However, that does not make sense when you get to the for loop. Your suggestion of rhs clearly works no matter how big xRange because I can get results for it even when the program falls over which just leaves me the question of why U1 = gather(AA)\gather(rhs); works and works quickly whereas U1 = AA\rhs; does not anbd gives NaN. Just before I started wrting this I thought that maybe I should look up and see what gather actually does.

gather - Transfer distributed array, Composite array or gpuArray to local workspace. This MATLAB function can operate on the following array data: On a gpuArray: transfers the elements of A from the GPU to the local workspace and assigns them to X.

I thought, because I skim read (it's an age thing!) that gather transferred stuff from GPU to ordinary memory. I'm not convinced that it actually does this. Couple of reasons. Firstly U1 is returned as a GPU array. If AA, a very large matrix, were being moved around I would expect to see something especially if gather is moving it to CPU RAM and then it would have to move the resultant back to put it into U1. I cannot see this in the diagrams previously shown. I just wondered if gather was a mechanism to provide a vector to the locations of the distributed data. This would be faster than looking for each piece individually. What is clear to me from the previous diagrams is that the CPU is NOT doing this calculation.

D = app.dt/(2*app.h^2); % D

A = diag(-2*ones(1,M)) + diag(ones(1,M-1),1) + diag(ones(1,M-1),-1);

A = sparse(A);

AD = sparse(D*A);

AS = sparse(1i*eye(size(A))+AD);%in this order I get the largest ...

% effective xRange and it is this line where finally OUT OF MEMORY ...

%causes a fail with size.

AB = gpuArray(complex(AD));

AA = gpuArray(complex(AS));

sdg = gpuArray(s*app.dt*app.gamma);%this is just a constant

clearvars A AD AS;

U = gpuArray(complex(zeros(M+2,app.NJ+1)));

U0 = gpuArray(complex(app.u0(2:end-1,:)));

for j=1:app.NJ

% This is the first shot at Crank Nicholson

rhs = (1i.*U0 - AB*U0 - sdg*(U0.*conj(U0)).*U0);

U1 = gather(AA)\gather(rhs); %r=1;

%which becomes a predictor corrector

U1 = (U1+U0)/2;

rhs = (1i.*U0 - AB*U0 - sdg*(U1.*conj(U1).*U1));

U1 = gather(AA)\gather(rhs); %r=2;

U0=U1;

U(:,j+1) = [0;U1;0];

end

[Bruno Luong]

"I thought, because I skim read (it's an age thing!) that gather transferred stuff from GPU to ordinary memory. I'm not convinced that it actually does this. Couple of reasons. Firstly U1 is returned as a GPU array."

That statement is incorrect, I set the break point right after

U1=gather(AA)\gather(rhs);

then do this

K>> whos U1
Name          Size             Bytes  Class     Attributes

U1        19999x1             319984  double    complex   

U1 is no longer gpuArray but standard double array, as you can see.

It is also what the gather doc says

"On a gpuArray: transfers the elements of A from the GPU to the local workspace and assigns them to X."

"it would have to move the resultant back to put it into U1"

No you missunderstand, U1 is on CPU after that. It is now assigned to a new DOUBLE array value. The old gpuArray class is irrelevant

[Joss Knight]

Do make sure you move the result of backslash back to the GPU so that subsequent operations take place on the GPU. So use gpuArray(gather(AA)\gather(rhs)) rather than just gather(AA)\gather(rhs), to ensure U1 is still a gpuArray on output.

As a developer of the gpuArray datatype, I'm usually - not always, but usually! - pretty confident in my assertions about its behaviour. Whatever you're seeing in the task manager is heavily confused by the fact that MATLAB's GPU functionality pools memory and operates lazily and asynchronously. That means the GPU isn't always being used when a line of code is run, it doesn't always finish when the next line of code is run, and memory no longer assigned to a variable is not necessarily released back to the system to show up as free in the Task Manager. Hopefully we will soon have our own monitor apps so you can more easily see this behaviour in MATLAB.

Worth reminding you perhaps that gather(AA) does not move AA to the CPU, it creates a temporary variable on the CPU and copies AA there. So the GPU memory is not released.

I'm glad you've found some compromises about where to do CPU work and where GPU. It takes quite a lot of memory to move a sparse CPU array to the GPU because a conversion between two different storage formats is required (CSC to CSR, if you care). It's usually best to create the data on the GPU from the outset, for instance by passing gpuArray data to spdiags.

[Bruno Luong]

Again, , what is the point of using gpuArray at all then, if the main calculation

gather(AA)\gather(rhs)

is done on cpu?

[Joss Knight]

The point is that every other calculation goes faster. Your numbers showed a big difference between CPU and GPU.

[Bruno Luong]

@Joss Knight But the major party runtime is when solving AA\rhs, not the rest. The advantage of solving linear system on GPU is lost when you force to do it on CPU.

The only reason for @Jonathan Wharrier to do CPU inversion because GPU is buggy beyond the spatial resolution of 54.

Here is the three methods run together, and clearly mixing between GPU and CPU (third result time) is the slowest, at least on my PC.

>> CPUtest
Config = cpuarray + home_native_backslash
	Elapsed time is 2.025 seconds
Config = gpuarray + home_native_backslash
	Elapsed time is 0.906 seconds
Config = gpuarray + cpu_backslash
	Elapsed time is 3.191 seconds

So still the question, what is the point? I mean you better off doing the whole calculation on CPU. The hybrid just slow down by 1.1 second, more than 50%.

The code demo is here.

close all;
configs = struct('arraytype',  {@cpuarray,             @gpuarray,             @gpuarray},...
                 'solvertype', {@home_native_backslash, @home_native_backslash, @cpu_backslash});
for k = 1:length(configs)
    array  = configs(k).arraytype;
    solver = configs(k).solvertype;
    % CRANK NICHOLSON VERSION
    xRange = 100;
    dx = 0.01; %mesh grid increment size
    x = -xRange:dx:xRange;%mesh grid
    T = 5; %time for run
    dt = 0.01; %Time increment size
    NJ=T/dt; %number of iterations
    t= dt*(0:NJ); %time vector
    M=length(x);    % length of the mesh
    M=M-2;          % active length minus end points
    S  = @(ex) 1*sech(sqrt(1/2)*ex).*exp(1i*2*ex); %creates Soliton
    u0 = S(x)'; %Start data - note transform
    s=1; %constant
    ul = zeros(size(t)); %set boundary values
    ur = zeros(size(t)); %boundaries are zero
    reset(gpuDevice);
    Bdy=zeros(M,NJ+1);
    A = spdiags(ones(3*M,1)*[1 -2 1],[-1 0 1],M,M);
    U = array(complex(zeros(M+2,NJ+1)));
    U0 = array(complex(u0(2:end-1,:)));
    UC = array(complex(size(U0)));
    D = dt/(2*dx^2); % D
    Bdy(1,:)= D*ul;
    Bdy(end,:)=D*ur;
    %Bdy = sparse(Bdy);
    url = array(complex(ul));
    urr = array(complex((ur)));
    Bdy = array(complex(Bdy));
    AB = array(complex(D*A));
    AA = 1i*speye(size(A))+AB;
    AA = array(AA);
    sdg = array(1*dt*1);
    t0=tic;
    for j=1:NJ
        CrankNicolsonRhsFun = @(U) (1i.*U0 - AB*U0 - sdg*(U.*conj(U)).*U - Bdy(:,j)-Bdy(:,j+1));
        rhs = CrankNicolsonRhsFun(U0);
        U1 = solver(AA,rhs);
        UC = (U1+U0)/2;
        rhs = CrankNicolsonRhsFun(UC);
        U1 = solver(AA,rhs);
        U0=U1;
        U(:,j+1) = [url(j+1);U1;urr(j+1)];
    end
    t=toc(t0);
    fprintf('Config = %s + %s\n', func2str(array), func2str(solver));
    fprintf('\tElapsed time is %1.3f seconds\n', t);
    U = gather(U);
    U(:,1) = double(u0);
    figure('Name', sprintf('Config(%d)', k));
    surf(abs(U),'EdgeColor','none');
end

function C = home_native_backslash(A,B)
C=A\B;
end

function C = cpu_backslash(A,B)
C=gpuArray(gather(A)\gather(B));
end

function a = cpuarray(a) % a input is double array
end

function a = gpuarray(a)
a = gpuArray(a);
end

[Jonathan Wharrier]

To make a run with xRange of 360 which I have done this morning means the process runs on a single thread. It is possible to use a multi threaded parallel approach, which I did in the early versions of the app. However there is a trade off. More threads, shared memory and therefore smaller xRange. The latter is unhelpful so one thread it is. If I run this solely on the CPU and watch task managerI found that I get a significant hump where the processor loads up the memory and performs the "for" calculation. Specifically on xR 300+ 4/min which is two passes through the Crank Nicholson/Predictor corrector calculation. I need a minimum of T = 25 with dt = 0.01 so 2500 goes through the CN/PC mill which takes approximately 20 hours. That is one run. I have a newer computer - still an AMD R7 but a 128Gb 5000 model rather than the 64GB R7 2400 laptop. It takes about an hour less on CPU.

This morning I have run so far, 10 runs. Each complete run T = 36 dt = 0.01, takes fractionally under 5 minutes. Now I do not claim to understand the why and the how but the difference on my machine is remarkable. Because I am pushing the machine to its limits - the CPU does not like 100% memory usage much and has fallen over twice - but I have a work around for the data - and the GPU falls over if it does not completely clear with "reset(gpuDevice)" and free up all the memory I have taken but again this seems to be an infrequent event, it is a bit frustrating BUT the overall is an improvement in the amount of data I am able to process by an order of magnitude.

[Joss Knight]

Agreed @Bruno Luong, this is what I see too, but it doesn't match Jonathan's figures which seemed to show an order of magnitude worse performance on CPU. If his figures had similar proportions to yours (and mine) I certainly wouldn't see much point in using the GPU, but also the original problem - that it's computationally intractable to run his code with high resolutions - is not there, since the CPU is only 2x slower. Most people are prepared to let things run twice as long if the alternative is something that occasionally crashes.

The GPU documentation is clear, it cannot always be faster, not all algorithms are straightforwardly parallelizable. If the matrix weren't tridiagonal, for instance, the GPU would be much slower than CPU. SVD is another algorithm that is highly problematic for GPU parallelization. Sparse in general is a peculiar problem, because GPU parallelization depends on uniformity but sparse data by its nature is non-uniform (except where it fits certain patterns like diagonality). That is why we don't provide a huge amount of sparse GPU functionality, because if it cannot be made faster than CPU we don't want people to use it only to be disappointed.

[Bruno Luong]

A "trick" of reduce cpu time is use decomposition object of AA matrix. Timing is now comparable to gpu.

Unfortunately there is no gpu support for decomposition.

>> CPUtest
Config = cpuarray + home_native_backslash
	Elapsed time is 1.980 seconds
Config = gpuarray + home_native_backslash
	Elapsed time is 0.896 seconds
Config = cpuarray + decomp
	Elapsed time is 1.156 seconds

with this demo script

close all;
configs = struct('arraytype',  {@cpuarray,              @gpuarray,              @cpuarray},...
                 'solvertype', {@home_native_backslash, @home_native_backslash, @decomp});
for k = 1:length(configs)
    array  = configs(k).arraytype;
    solver = configs(k).solvertype;
    % CRANK NICHOLSON VERSION
    xRange = 100;
    dx = 0.01; %mesh grid increment size
    x = -xRange:dx:xRange;%mesh grid
    T = 5; %time for run
    dt = 0.01; %Time increment size
    NJ=T/dt; %number of iterations
    t= dt*(0:NJ); %time vector
    M=length(x);    % length of the mesh
    M=M-2;          % active length minus end points
    S  = @(ex) 1*sech(sqrt(1/2)*ex).*exp(1i*2*ex); %creates Soliton
    u0 = S(x)'; %Start data - note transform
    s=1; %constant
    ul = zeros(size(t)); %set boundary values
    ur = zeros(size(t)); %boundaries are zero
    reset(gpuDevice);
    Bdy=zeros(M,NJ+1);
    A = spdiags(ones(3*M,1)*[1 -2 1],[-1 0 1],M,M);
    U = array(complex(zeros(M+2,NJ+1)));
    U0 = array(complex(u0(2:end-1,:)));
    UC = array(complex(size(U0)));
    D = dt/(2*dx^2); % D
    Bdy(1,:)= D*ul;
    Bdy(end,:)=D*ur;
    %Bdy = sparse(Bdy);
    url = array(complex(ul));
    urr = array(complex((ur)));
    Bdy = array(complex(Bdy));
    AB = array(complex(D*A));
    AA = 1i*speye(size(A))+AB;
    AA = array(AA);
    sdg = array(1*dt*1);
    t0=tic;
    if isequal(solver, @decomp)
        AA= decomposition(AA);
    end
    for j=1:NJ
        CrankNicolsonRhsFun = @(U) (1i.*U0 - AB*U0 - sdg*(U.*conj(U)).*U - Bdy(:,j)-Bdy(:,j+1));
        rhs = CrankNicolsonRhsFun(U0);
        U1 = solver(AA,rhs);
        UC = (U1+U0)/2;
        rhs = CrankNicolsonRhsFun(UC);
        U1 = solver(AA,rhs);
        U0=U1;
        U(:,j+1) = [url(j+1);U1;urr(j+1)];
    end
    t=toc(t0);
    fprintf('Config = %s + %s\n', func2str(array), func2str(solver));
    fprintf('\tElapsed time is %1.3f seconds\n', t);
    U = gather(U);
    U(:,1) = double(u0);
    figure('Name', sprintf('Config(%d)', k));
    surf(abs(U),'EdgeColor','none');
end

function C = home_native_backslash(A,B)
C=A\B;
end

function C = decomp(A,B)
C=A\B; % A is decompistion class, not matrix
end

function C = cpu_backslash(A,B)
C=gpuArray(gather(A)\gather(B)); % what is the point compare to purely cpu?
end

function a = cpuarray(a) % a input is double array
end

function a = gpuarray(a)
a = gpuArray(a);
end

[Jonathan Wharrier]

Words fail me frankly. I am deeply impressed..At the moment I need to do some writing up and, because of fast gpu runs I now have enough runs in the short term. However, you have made me think and I will file this away for "play" later this year. Like you I am puzzled by C=gpuArray(gather(A)\gather(B)); I haven't put the gpuArray in yet and nothing on taskManager shows CPU activity during the CN solve! I also went back to my old code using GPU just to check and it still is 100+ times slower.

The long and the short is that I have now a sensible operable speed for the calculations and have learned a great deal for which I am very grateful to yourself (and also Joss.).

[Bruno Luong]

@Jonathan Wharrier May be if you are satisfied with the answer and it somehow sheds a light even if it doesn't still explain few obscure aspects, could you accept the answer?

[Jonathan Wharrier]

I am happy to do so. I have found the whole experience very helpful.

Tx Jo