[GSOC24] Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2 #2346

areenraj · 2024-08-27T10:01:24Z

Proposed Changes

This is the modified version of SU2 code that supports CUDA usage for the FGMRES solver and the use of NVBLAS. The main focus is the offloading of the Matrix Vector Product in the FGMRES solver to the GPU using CUDA Kernels. This implementation shows promise with marginally better run times (all benchmarks were carried out with the GPU Error Checking switched off and in debug mode to check if the correct functions were being called).

The use of NVBLAS is secondary and while functionality has been added to make it usable, it is not activated as it doesn't cause an appreciative increase in performance.

Compilation and Usage

Compile using the following MESON Flag

-Denable-cuda=true

And activate the functions using the following Config File Option

ENABLE_CUDA=YES

NOTE ON IMPLEMENTATION

I've decided to go with a single version of the code where the CPU and GPU implementations co-exist in the same linear solver and can be disabled or switched using a combination of Meson and Config File options. This is why I have defined three classes - one over-arching class that is named CExecutionPath that has two child classes - CCpuExecution and CGpuExecution. These child classes contain the correct member function for each path - CPU or GPU functioning.

All of this could also be easily achieved with an if statement that switches between the two - but that particular implementation will access and run the statement for each call. In our case once a Matrix Vector Product object is created, it will immediately know whether to use CPU or GPU mode of execution.

Recommendations are most welcome to improve or make this implementation better

PR Checklist

Warning Levels do come (only at level 3) but they are all of the following type

style of line directive is a GCC extension

The documentation for compiling with CUDA needs to be added by forking the SU2 site repo and adding the relevant changes to it? Or do I need to contact someone to change things on the site itself?

DOxygen Documentations and config template are all updated.

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

…on as well

Created Final Report

pcarruscag

Nice work

pcarruscag · 2024-08-27T12:04:27Z

Common/src/linear_algebra/GPU_lin_alg.cu

+
+  gpuErrChk(cudaMemcpy((void*)(d_matrix), (void*)&matrix[0], (sizeof(ScalarType)*mat_size), cudaMemcpyHostToDevice));
+  gpuErrChk(cudaMemcpy((void*)(d_vec), (void*)&vec[0], (sizeof(ScalarType)*vec_size), cudaMemcpyHostToDevice));
+  gpuErrChk(cudaMemcpy((void*)(d_prod), (void*)&prod[0], (sizeof(ScalarType)*vec_size), cudaMemcpyHostToDevice));


you don't need to copy the product, you just need to memset to 0

Good catch, will add this. Thank you

pcarruscag · 2024-08-27T12:05:20Z

Common/src/linear_algebra/GPU_lin_alg.cu

+  double xDim = (double) 1024.0/(nVar*nEqn);
+  dim3 blockDim(floor(xDim), nVar, nEqn);
+  double gridx = (double) nPointDomain/xDim;
+  dim3 gridDim(ceil(gridx), 1, 1);
+


Can you document the choice of work distribution between blocks and threads?

Full report with Benchmark Results and other Explanations

Complete definition of how the Kernel is initialized

Explanation as to why the blocks are multi-directional

Algorithm for Block Matrix Calculations

pcarruscag · 2024-08-27T12:08:29Z

Common/src/linear_algebra/GPU_lin_alg.cu

+      for(int index = d_row_ptr[i]; index<d_row_ptr[i+1]; index++)
+      {
+        int matrix_index = index * nVar * nEqn;
+        int vec_index = d_col_ind[index] * nEqn;
+
+        res += matrix[matrix_index + (j * nEqn + k)] * vec[vec_index + k];
+      }


Is this based on some publication? Did you experiment with other divisions of work?
For example, I see you are going for coalesced access to the matrix blocks, but this requires multiple reads of the same vector entries.

I haven't experimented with different implementations as I went with this because it seemed optimal. It does access the same vector elements repeatedly while going through an entire row.

Haven't particularly checked out publications yet but you're right, there may be a better way to do this and I'll look into them. If you do have some recommendations on improvements then let me know.

Also, the current bottleneck is based on memory copy between the CPU and GPU while the kernel launch itself is 20 times faster than the repeated copy. Any insights on that as well would be very grateful.

Our current approach to circumvent this is to port not only single subroutines like matrix vector multiplication but the entire Krylov Solver loop where it searches over all the search directions in the subspace to the GPU. This would cut down on the repeated memory transfers.

You can do something like a flag in CSysMatrix that is set to true when the matrix is uploaded to GPU and set to false when the matrix changes (for example when we clear the matrix to write new blocks).
You can use pinned host memory to make the transfers faster.
You can try uploading the matrix in chunks and overlap the uploads with the CPU work of filling another chunk.
Ultimately, the issue of transferring the matrix only goes away by porting the entire code 😅

Regarding recommendations I have to be a bit cryptic because of my current job, but the general goals are coalesced access, and avoid reading or writing the same global memory location more than once.
I read this paper before my current job.
Optimization of Block Sparse Matrix-Vector Multiplication on Shared-Memory
Parallel Architectures
But like you said, there is much more to gain by porting more linear algebra operations than to micro-optimize the multiplications.

Hmm, I'm a bit confused here. My entire project and the previous pull request's purpose has been to slowly transfer the FGMRES solver to GPU.

Do you mean that we should just port the helper functions like Matrix Vector Product and other linear algebra functions and let the solver decide based on user given parameters?

Also I'll try to catch a dev meeting so that we can maybe discuss how to move forward from here.

Correct, the same way the FGMRES doesn't need to know how to multiply matrices and vectors, it doesn't need to know anything about GPUs

Got it, I will get on this and will also discuss some things I have in mind during the dev meet

I can't attend the dev meeting today, I suggest you open an issue and that way we can talk about the plan and the GSOC candidate that may work on the project also benefit from us writing everything down.

...This should increase performance (although I am still testing it).

NSight shows the new Kernel is twice as fast. Although this doesn't mean anything for overall performance as the bottleneck is still memory transfer.

kursatyurt · 2024-10-07T18:24:04Z

I could not follow the work here; I just watched the conference presentation. As much work has been done, have you considered using an existing linear algebra library such as https://github.com/ginkgo-project/ginkgo?

AFAIK the only thing you need to do is copy the matrix to Gingko format on GPU. Then Gingko will provide an efficient, scalable solver that works not only on NVIDIA but also on AMD and Intel.

GMRES and ILU preconditioners are available there, so it is pretty much ready to go for all problems.

areenraj · 2024-10-07T18:46:09Z

@kursatyurt Hello, thank you so much for the lead.

Our initial scope mostly involved writing our own kernels and I did explore some libraries at the start - I was planning on using CUSP as well but my main concern was its lack of being updated to the newly compatible versions of the toolkit. cuSolver and cuBLAS do exist, but I chose to go ahead with a "simple" kernel implementation to have more control. I also felt that if I could keep the block size of the grid in optimal territory then they could be just as fast as those options (please do correct me if my reading of the literature or the situation was incorrect)

I was not aware of Ginkgo and I will surely give it a go and try to produce some comparative results. I am currently super busy for this month and will get to working on the code with some delay.

Again, thank you for the lead!

areenraj · 2024-10-07T18:47:06Z

Also could you mention what you meant by "you could not follow the work here" If there is a specific doubt in the work then I would love to clarify it over slack whenever I get the time

kursatyurt · 2024-10-07T19:41:46Z

Also could you mention what you meant by "you could not follow the work here" If there is a specific doubt in the work then I would love to clarify it over slack whenever I get the time

I am not very familiar with the linear solver implementation in SU2 not about the work itself

kursatyurt · 2024-10-07T20:12:53Z

@kursatyurt Hello, thank you so much for the lead.

Our initial scope mostly involved writing our own kernels and I did explore some libraries at the start - I was planning on using CUSP as well but my main concern was its lack of being updated to the newly compatible versions of the toolkit. cuSolver and cuBLAS do exist, but I chose to go ahead with a "simple" kernel implementation to have more control. I also felt that if I could keep the block size of the grid in optimal territory then they could be just as fast as those options (please do correct me if my reading of the literature or the situation was incorrect)

To learn the basics, it's a good idea, but for large-scale projects, I prefer using existing libraries if possible.
Those libraries generally exploit state-of-the-art solution like mixed-precision computing. A gaming GPU is not way faster than a good CPU in double precision, but way faster in single precision, most of them have 64:1 ratio, however server class GPU have 2:1 ratio. Also when available they use vendor libraries like cuBLAS or hipBLAS. It is always nice to have you only care about connection and somebody else handle the solver as performant as possible. In future probably they will provide more and more solvers and it will be automagically works.

It is kind of light-weight too, not a huge dependency like Trilinos or PETSc.

I was not aware of Ginkgo and I will surely give it a go and try to produce some comparative results. I am currently super busy for this month and will get to working on the code with some delay.

Again, thank you for the lead!

I can test on various GPUs (P100/V100/A100 and 4070Mobile) on single node multi-gpu etc.

pcarruscag

Is your idea to return to this project as part of GSOC or just personal interest?

pcarruscag · 2025-03-17T13:03:06Z

Common/src/linear_algebra/CSysMatrix.cpp

+#if defined(HAVE_CUDA)
+  gpuErrChk(cudaMalloc((void**)(&d_row_ptr), (sizeof(row_ptr) * (nPointDomain + 1.0))));
+  gpuErrChk(cudaMalloc((void**)(&d_col_ind), (sizeof(col_ind) * nnz)));
+  gpuErrChk(cudaMallocHost((void**)(&matrix), (sizeof(ScalarType) * nnz * nVar * nEqn)));


what is the idea with this change?

This prevents cuda Memory functions like Malloc to be called by the compiler when the user doesn't have CUDA. The conditional is only executed when the user specifies that they would like to compile with cuda and have the toolkit along with nvcc installed. If I leave it as is then when compiling with just gcc will throw an exception.

If there is a different/better way of doing this then please do let me know :D

I'm talking about cudaMallocHost instead of cudaMalloc

Oh that was done to store the matrix in pinned memory since transfers are faster to it

But you removed d_matrix, so now the GPU is working directly from host memory which is probably slow.
Also if you read the docs for cudaMallocHost it's not advised to use large chunks of pinned memory, and the matrix is the largest block of memory in SU2...

But you removed d_matrix, so now the GPU is working directly from host memory which is probably slow.

My reasoning was that when pinnned memory is allocated on 64 bit OSes, the memory is also simultaneously mapped to the GPU as the conditions for universal virtual addressing are met. This method proves to be faster than the regular MemCpy when the data address is only read/written to once. In our case, the Mat Vec Kernel only read the matrix values once per CudaMemCpy - so I switched it out this way. If my understanding of this was wrong then do let me know.

But now that we will be working on porting the other linear algebra functions it will make more sense to use the CudaMemCpy approach of transferring the pinned memory to the GPU. As the same variables will be accessed multiple times over the course of a single memory transfer

Also if you read the docs for cudaMallocHost it's not advised to use large chunks of pinned memory, and the matrix is the largest block of memory in SU2

i never ran into a problem with this as the cases I ran didn't take a performance hit. But I can understand that with increase in problem size the pinned memory could decrease the total available system memory and bottleneck the performance of the solver. Would you recommend using pinned memory for the input and output vectors then?

areenraj · 2025-03-17T15:45:50Z

Is your idea to return to this project as part of GSOC or just personal interest?

@pcarruscag No, I won't be applying to GSoC any time in the future. Its just that I started a work and I would like to see this through. If anyone is selected for GPU Acceleration this time then I would be more than happy to work and collaborate with them to make sure the project receives as much help as it can.

…ansfer and standardized file structure

areenraj · 2025-03-29T17:27:51Z

29 Mar 2025 Commits

The commits currently made aim to provide more control over memory access and transfer so that we can make further changes to the surrounding linear algebra functions which make up the FGMRES Solver.

Streamlined Vector Storage in the GPU

Each vector definition now also creates its corresponding analog in the GPU Memory space. This is done by allocating memory for it using cudaMalloc when the Initialize function of the CSysVector Class is called. This allows for the continuous storage of vector data in the GPU memory in between calls of the Matrix Vector Product and other linear algebra functions. The previous implementation only allowed this data to persist during a single call of the Matrix Vector Product and it had to be refreshed for each call.

This implementation was similar to how the matrix was stored previously. Saving the matrix in pinned memory is also removed due to its huge size as pointed out by earlier feedback. The device pointer can be accessed at any point using a new dedicated public member function (GetDevicePointer) of the CsysVector Class.

Added Memory Transfer Control

As previously discussed, we needed more control over the memory transfer in between calls so that a larger load of the computation could be carried out on the GPU without multiple memory transfers. Now these transfers are carried out by member functions with a flag built into them to decide whether the copy needs to be carried out or not. This flag is set to true by default and does not need to be specified all the time.

Further changes are necessary to actually use this flag to decrease the frequency of memory transfer - namely a variable that allows the inner loop of FGMRES to communicate with the MatrixVectorProduct function to know when to switch the flag on or off. This will be added after I port the preconditioner.

Minor Change - Standardized File Structure Slightly

Redundant .cuh header files are now gone. I have added a GPUComms.cuh file so that any functions that need to be accessed by all CUDA files - like error checking - can be added here for future references. I've also added GPUVector and GPUMatrix files here - each containing the cuda wrapping member functions for the CSysVector and CSysMatrix class respectively.

Please let me know if you notice any bugs or if my implementations can be improved in any way.

areenraj · 2025-03-29T17:28:36Z

I'll now move forward and start porting the LU_SGS Preconditioner

pcarruscag · 2025-03-30T16:07:05Z

Common/src/linear_algebra/GPUVector.cu

Rename to CSysVectorGPU or something similar, this implements functions of CSysVector so it makes more sense to name it consistently with the CPU version.

pcarruscag

Can you fix the merge conflicts so the regression tests can run in this PR?

pcarruscag · 2025-03-30T16:10:03Z

Common/include/linear_algebra/CMatrixVectorProduct.hpp

+    std::cerr << "\nError in launching Matrix-Vector Product Function\n";
+    std::cerr << "ENABLE_CUDA is set to YES\n";
+    std::cerr << "Please compile with CUDA options enabled in Meson to access GPU Functions" << std::endl;


SU2_MPI::Error(...

This is fixed, however in the gpuErrChk function I have left it as fprintf since the the string has some placeholders to it. I could use a temporary character array along with SPRINTF to use SU2_MPI::Error. But I am not sure whether allocating and de-allocating the temporary character array repeatedly will slow down performance or not. Since, gpuErrChk is called every time a CUDA API function is called.

pcarruscag · 2025-03-30T16:11:23Z

Common/include/linear_algebra/CSysMatrix.hpp

+
+  void FGMRESMainLoop(std::vector<ScalarType> W, std::vector<ScalarType> Z, su2vector<ScalarType>& g,
+                      su2vector<ScalarType>& sn, CSysVector<ScalarType>& cs, su2vector<ScalarType>& y,
+                      su2vector<ScalarType>& H, int m, CGeometry* geometry, const CConfig* config) const;


This seems unused

Suggested change

void FGMRESMainLoop(std::vector<ScalarType> W, std::vector<ScalarType> Z, su2vector<ScalarType>& g,

su2vector<ScalarType>& sn, CSysVector<ScalarType>& cs, su2vector<ScalarType>& y,

su2vector<ScalarType>& H, int m, CGeometry* geometry, const CConfig* config) const;

pcarruscag · 2025-03-30T16:14:36Z

Common/src/linear_algebra/CSysVector.cpp

  if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc<ScalarType, true>(64, nElm * sizeof(ScalarType));

+#if defined(HAVE_CUDA)
+  gpuErrChk(cudaMalloc((void**)(&d_vec_val), (sizeof(ScalarType) * nElm)));
+  gpuErrChk(cudaMemset((void*)d_vec_val, 0.0, (sizeof(ScalarType) * nElm)));


Similarly to how we have MemoryAllocation::aligned_alloc and free, please create wrappers for the cuda-specific functions so that we avoid #ifdef HAVE_CUDA in many places.

Done, please let me know if my handling of null pointers was accurate

pcarruscag · 2025-03-30T16:15:07Z

Common/src/linear_algebra/GPUMatrix.cu

+   int row = (blockIdx.x * blockDim.x + threadIdx.x)/32;
+   int threadNo = threadIdx.x%32;


Define 32 (threads per warp) as a constant somewhere.

pcarruscag · 2025-03-30T16:16:27Z

Common/src/linear_algebra/GPUMatrix.cu

+  dim3 blockDim(1024,1,1);
+  double gridx = (double) nPointDomain/32.0;
+  gridx = double(ceil(gridx));


We have functions to round up integer division in the memory allocation toolbox.
Define 1024 as the default block size somewhere too.

Used a slight variation of the function in the Allocation Toolbox. Similar syntax however.

update master

…k Sizes

areenraj · 2025-05-08T15:32:11Z

Can you fix the merge conflicts so the regression tests can run in this PR?

It shows that merging is blocked (I guess a review is required?) Also please let me know how you wish for me to keep the fork synced with changes in the original - either through merge or rebase.

areenraj added 30 commits July 26, 2024 17:53

refresh everything

34c0bda

readme update

c953af8

final push

231968d

Enable GPU Mat Vec

4bddc30

New Branch and Optimized Memory Alloc on GPU Slightly

2da57f6

Finished GPU Mat-Vec with CPU Accuracy and Block Matrix Parallelizati…

aa35779

…on as well

Added fully template kernels to prevent type errors

8866936

Disabled NVBLAS Implementation in the DG Solver

12f01db

Added options and error check

d7cbd5e

reverting stuff to debug

103ec39

fixed turbulent case but Error Check is a performance hit

38d658b

Updated README for final report

02b9eb8

readme update

1a902ae

readme update

3fa00a9

readme update

57ffb74

Final Graph Changes

b8f14d3

Image Changes

3a000e2

Added Runtime Polymorphism for selecting execution Path

c31b9df

Added Runtime Polymorphism to select between CPU and GPU Execution

e131308

Added runtime polymorphism to select execution path

7695314

Added runtime polymorphism to select execution path

809c2d0

Added Runtime Polymorphism to select between CPU and GPU Execution

9832f33

Added Runtime Polymorphism to select between CPU and GPU Execution

1f41592

Create REPORT.md

de0810a

Created Final Report

Delete REPORT.md

66878d1

Added Preprocessor Directives

729bfc8

Making Repo PR Ready

3f351db

Making Repo PR Ready

f363809

Making it PR Ready

a0e09d7

Pre-Commit Hook Ran

b489d09

PR Ready

4926d34

areenraj changed the title ~~Addition of CUDA and GPU Execution Capabilities to FGMRES Linear Solver in SU2~~ Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2 Aug 27, 2024

pcarruscag reviewed Aug 27, 2024

View reviewed changes

areenraj changed the title ~~Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2~~ [GSOC24] Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2 Aug 29, 2024

added some fixes and error handling

c691960

areenraj added 3 commits March 17, 2025 13:32

New Algorithm for MVP

1cc128d

Merge remote-tracking branch 'upstream/master' into HEAD

9fb99b2

Version Number Change

86e59b4

pcarruscag reviewed Mar 17, 2025

View reviewed changes

Streamlined Vector GPU Storage, Introduced flags to control memory tr…

c7777c8

…ansfer and standardized file structure

pcarruscag reviewed Mar 30, 2025

View reviewed changes

pcarruscag and others added 4 commits May 4, 2025 16:27

Merge pull request su2code#2499 from su2code/develop

5d71ae2

update master

Introduced Memory Wrapper Functions, Error Handling and Constant Bloc…

b7ed530

…k Sizes

Syncing Repo using Merge

625c483

Minor changes

70a7680

		int row = (blockIdx.x * blockDim.x + threadIdx.x)/32;
		int threadNo = threadIdx.x%32;

[GSOC24] Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2 #2346

Are you sure you want to change the base?

[GSOC24] Addition of CUDA and GPU Acceleration to FGMRES Linear Solver in SU2 #2346

Uh oh!

Conversation

areenraj commented Aug 27, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Compilation and Usage

NOTE ON IMPLEMENTATION

PR Checklist

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj Aug 29, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj Aug 29, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj Aug 29, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj Mar 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kursatyurt commented Oct 7, 2024

Uh oh!

areenraj commented Oct 7, 2024

Uh oh!

areenraj commented Oct 7, 2024

Uh oh!

kursatyurt commented Oct 7, 2024

Uh oh!

kursatyurt commented Oct 7, 2024

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj Mar 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

areenraj commented Mar 17, 2025

Uh oh!

areenraj commented Aug 27, 2024 •

edited

Loading

areenraj Aug 29, 2024 •

edited

Loading

areenraj Aug 29, 2024 •

edited

Loading

areenraj Aug 29, 2024 •

edited

Loading

areenraj Mar 17, 2025 •

edited

Loading

areenraj Mar 17, 2025 •

edited

Loading