From 2fdf7aa8c1ff8a651fe385b0a6672782b14501f2 Mon Sep 17 00:00:00 2001 From: Matthew Nicely Date: Tue, 11 May 2021 17:13:26 -0400 Subject: [PATCH 1/3] added streams and cufft --- 2D_Program/Makefile | 2 +- 2D_Program/solver_functions/DST.c | 103 +++++++++++------ 2D_Program/solver_functions/cuda_kernels.cu | 74 ++++++------ 2D_Program/solver_functions/cuda_kernels.h | 26 ++++- 2D_Program/solver_functions/solver.c | 121 ++++++++++---------- 5 files changed, 194 insertions(+), 132 deletions(-) diff --git a/2D_Program/Makefile b/2D_Program/Makefile index 1112719..bfb90f7 100644 --- a/2D_Program/Makefile +++ b/2D_Program/Makefile @@ -28,7 +28,7 @@ endif ifneq ($(origin USE_CUFFTW), undefined) ifneq ($(USE_CUFFTW), 0) CFLAGS += -DUSE_CUFFTW - LIBS += -lcufft -lcufftw -lcudart + LIBS += -lcufft -lcudart else LIBS += -lfftw3_omp -lfftw3 endif diff --git a/2D_Program/solver_functions/DST.c b/2D_Program/solver_functions/DST.c index 8b95c5d..c6cb2f2 100644 --- a/2D_Program/solver_functions/DST.c +++ b/2D_Program/solver_functions/DST.c @@ -8,44 +8,81 @@ #if USE_CUFFTW #ifdef USE_COMBINE void fullDST( const cudaStream_t *streams, + const cudaEvent_t * events, const System sys, const DSTN dst, - const fftw_plan plan, - const fftw_plan plan2, + const cufftHandle plan, + cuDoubleComplex * d_workspace[2], cuDoubleComplex * d_y, double * in, fftw_complex * out, double * in2, fftw_complex * out2 ) { - PUSH_RANGE( "forwardDST", 2 ) - CUDA_RT_CALL( cudaStreamSynchronize( streams[0] ) ); - load_1st_DST_wrapper( sys, dst, sys.rhs, in, in2 ); + PUSH_RANGE( "1st DST", 2 ) + CUDA_RT_CALL( cudaEventRecord( events[2], streams[2] ) ); - CUDA_RT_CALL( cudaStreamSynchronize( streams[1] ) ); - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan creation + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 creation + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait for sys.rhs + + load_1st_DST_wrapper( streams[0], sys, dst, sys.rhs, in, in2 ); + + CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); + + CUDA_RT_CALL( cudaStreamWaitEvent( streams[1], events[0], cudaEventWaitDefault ) ); // Wait for load_1st_DST + CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); + CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[0] ) ); + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + + CUDA_RT_CALL( cufftSetStream( plan, streams[1] ) ); + CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[1] ) ); + CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] + + CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); + CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); + CUDA_RT_CALL( cudaEventRecord( events[2], streams[3] ) ); POP_RANGE - PUSH_RANGE( "forwardDST", 3 ) - CUDA_RT_CALL( cudaStreamSynchronize( streams[2] ) ); - middle_stuff_ls_DST_wrapper( sys, dst, out, out2, in, in2, d_y ); + PUSH_RANGE( "Trig Solver", 3 ) + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan execution + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 execution + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait forsys.U, sys.L, sys.Up + + middle_stuff_ls_DST_wrapper( streams[0], sys, dst, out, out2, in, in2, d_y ); + + CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); POP_RANGE - PUSH_RANGE( "forwardDST", 4 ) - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ + PUSH_RANGE( "2nd DST", 4 ) + CUDA_RT_CALL( cudaStreamWaitEvent( streams[1], events[0], cudaEventWaitDefault ) ); // Wait for middle_stuff_ls_DST + CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); - CUDA_RT_CALL( cudaStreamSynchronize( streams[3] ) ); - store_2st_DST_wrapper( sys, dst, out, out2, sys.sol ); + CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[0] ) ); + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + + CUDA_RT_CALL( cufftSetStream( plan, streams[1] ) ); + CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[1] ) ); + CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] + + CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); + CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); + CUDA_RT_CALL( cudaEventRecord( events[2], streams[4] ) ); + + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan execution + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 execution + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait for sys.sol + store_2st_DST_wrapper( streams[0], sys, dst, out, out2, sys.sol ); + + CUDA_RT_CALL( cudaStreamSynchronize( streams[0] ) ); // Wait for store_2st_DST POP_RANGE } #else void fullDST( const cudaStream_t *streams, + const cudaEvent_t * events, const System sys, const DSTN dst, - const fftw_plan plan, - const fftw_plan plan2, + const cufftHandle plan, cuDoubleComplex * d_rhat, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y, @@ -54,28 +91,26 @@ void fullDST( const cudaStream_t *streams, double * in2, fftw_complex * out2 ) { - PUSH_RANGE( "forwardDST", 2 ) - CUDA_RT_CALL( cudaStreamSynchronize( streams[0] ) ); - load_1st_DST_wrapper( sys, dst, sys.rhs, in, in2 ); + PUSH_RANGE( "1st DST", 2 ) + load_1st_DST_wrapper( NULL, sys, dst, sys.rhs, in, in2 ); + + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] - CUDA_RT_CALL( cudaStreamSynchronize( streams[1] ) ); - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ - store_1st_DST_wrapper( sys, dst, out, out2, d_rhat ); + store_1st_DST_wrapper( NULL, sys, dst, out, out2, d_rhat ); POP_RANGE - PUSH_RANGE( "forwardDST", 3 ) - CUDA_RT_CALL( cudaStreamSynchronize( streams[2] ) ); - middle_stuff_DST_wrapper( sys, d_rhat, d_xhat, d_y ); + PUSH_RANGE( "Trig Solver", 3 ) + middle_stuff_DST_wrapper( NULL, sys, d_rhat, d_xhat, d_y ); POP_RANGE - PUSH_RANGE( "forwardDST", 4 ) - CUDA_RT_CALL( cudaStreamSynchronize( streams[3] ) ); - load_2st_DST_wrapper( sys, dst, d_xhat, in, in2 ); - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ + PUSH_RANGE( "2nd DST", 4 ) + load_2st_DST_wrapper( NULL, sys, dst, d_xhat, in, in2 ); + + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] - store_2st_DST_wrapper( sys, dst, out, out2, sys.sol ); + store_2st_DST_wrapper( NULL, sys, dst, out, out2, sys.sol ); POP_RANGE } #endif diff --git a/2D_Program/solver_functions/cuda_kernels.cu b/2D_Program/solver_functions/cuda_kernels.cu index eb0d0c3..a4327ac 100644 --- a/2D_Program/solver_functions/cuda_kernels.cu +++ b/2D_Program/solver_functions/cuda_kernels.cu @@ -206,24 +206,24 @@ __global__ void __launch_bounds__( 256 ) middle_stuff_DST( const int N, // } // } __global__ void __launch_bounds__( 64 ) middle_stuff_ls_DST( const int N, - const int Nx, - const int Ny, - const int NC, - const double coef, - const cuDoubleComplex *__restrict__ out, - const cuDoubleComplex *__restrict__ out2, - const cuDoubleComplex *__restrict__ d_SysU, - const cuDoubleComplex *__restrict__ d_SysL, - const cuDoubleComplex *__restrict__ d_SysUp, - cuDoubleComplex *__restrict__ d_y, - double *__restrict__ in, - double *__restrict__ in2 ) { + const int Nx, + const int Ny, + const int NC, + const double coef, + const cuDoubleComplex *__restrict__ out, + const cuDoubleComplex *__restrict__ out2, + const cuDoubleComplex *__restrict__ d_SysU, + const cuDoubleComplex *__restrict__ d_SysL, + const cuDoubleComplex *__restrict__ d_SysUp, + cuDoubleComplex *__restrict__ d_y, + double *__restrict__ in, + double *__restrict__ in2 ) { const int tidX { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; cuDoubleComplex temp {}; - if (tidX < Nx) { + if ( tidX < Nx ) { int mx = Ny * tidX; temp = make_cuDoubleComplex( -out[tidX + 1].y, -out2[tidX + 1].y ); @@ -242,7 +242,7 @@ __global__ void __launch_bounds__( 64 ) middle_stuff_ls_DST( const int N, in[( Ny - 1 ) * N + tidX + 1] = temp.x; in2[( Ny - 1 ) * N + tidX + 1] = temp.y; -#pragma unroll 8 +#pragma unroll 4 for ( int j = Ny - 2; j >= 0; j-- ) { cuDoubleComplex temp2 = cuCdiv( cuCsub( d_y[j * Ny + tidX], cuCmul( d_SysUp[mx + j], temp ) ), d_SysU[mx + j] ); @@ -253,7 +253,12 @@ __global__ void __launch_bounds__( 64 ) middle_stuff_ls_DST( const int N, } } -void load_1st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleComplex *d_rhs, double *in, double *in2 ) { +void load_1st_DST_wrapper( const cudaStream_t stream, + const System sys, + const DSTN dst, + const cuDoubleComplex *d_rhs, + double * in, + double * in2 ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -266,13 +271,13 @@ void load_1st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleCompl void *args[] { &N, &Nx, &Ny, &d_rhs, &in, &in2 }; - CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } -void store_1st_DST_wrapper( const System sys, +void store_1st_DST_wrapper( const cudaStream_t stream, + const System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, @@ -291,13 +296,17 @@ void store_1st_DST_wrapper( const System sys, void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &d_rhat }; - CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } -void load_2st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleComplex *d_xhat, double *in, double *in2 ) { +void load_2st_DST_wrapper( const cudaStream_t stream, + const System sys, + const DSTN dst, + const cuDoubleComplex *d_xhat, + double * in, + double * in2 ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -310,13 +319,13 @@ void load_2st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleCompl void *args[] { &N, &Nx, &Ny, &d_xhat, &in, &in2 }; - CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } -void store_2st_DST_wrapper( const System sys, +void store_2st_DST_wrapper( const cudaStream_t stream, + const System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, @@ -335,13 +344,13 @@ void store_2st_DST_wrapper( const System sys, void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &d_sol }; - CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } -void middle_stuff_DST_wrapper( System sys, +void middle_stuff_DST_wrapper( const cudaStream_t stream, + System sys, const cuDoubleComplex *d_rhat, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y ) { @@ -357,13 +366,13 @@ void middle_stuff_DST_wrapper( System sys, void *args[] { &N, &Nx, &Ny, &sys.U, &sys.L, &sys.Up, &d_rhat, &d_xhat, &d_y }; - CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &middle_stuff_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &middle_stuff_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } -void middle_stuff_ls_DST_wrapper( System sys, +void middle_stuff_ls_DST_wrapper( const cudaStream_t stream, + System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, @@ -380,13 +389,12 @@ void middle_stuff_ls_DST_wrapper( System sys, CUDA_RT_CALL( cudaDeviceGetAttribute( &numSMs, cudaDevAttrMultiProcessorCount, 0 ) ); int threadPerBlock { 64 }; - int blocksPerGrid { (N + threadPerBlock - 1) / threadPerBlock }; + int blocksPerGrid { ( N + threadPerBlock - 1 ) / threadPerBlock }; void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &sys.U, &sys.L, &sys.Up, &d_y, &in, &in2 }; CUDA_RT_CALL( - cudaLaunchKernel( ( void * )( &middle_stuff_ls_DST ), blocksPerGrid, threadPerBlock, args, 0, NULL ) ); + cudaLaunchKernel( ( void * )( &middle_stuff_ls_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); CUDA_RT_CALL( cudaPeekAtLastError( ) ); - CUDA_RT_CALL( cudaStreamSynchronize( NULL ) ); } \ No newline at end of file diff --git a/2D_Program/solver_functions/cuda_kernels.h b/2D_Program/solver_functions/cuda_kernels.h index 759ce73..6715ebf 100644 --- a/2D_Program/solver_functions/cuda_kernels.h +++ b/2D_Program/solver_functions/cuda_kernels.h @@ -6,23 +6,37 @@ extern "C" { #endif -void load_1st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleComplex *d_rhs, double *in, double *in2 ); -void store_1st_DST_wrapper( const System sys, +void load_1st_DST_wrapper( const cudaStream_t streams, + const System sys, + const DSTN dst, + const cuDoubleComplex *d_rhs, + double * in, + double * in2 ); +void store_1st_DST_wrapper( const cudaStream_t streams, + const System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, cuDoubleComplex * d_rhat ); -void load_2st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleComplex *d_xhat, double *in, double *in2 ); -void store_2st_DST_wrapper( const System sys, +void load_2st_DST_wrapper( const cudaStream_t streams, + const System sys, + const DSTN dst, + const cuDoubleComplex *d_xhat, + double * in, + double * in2 ); +void store_2st_DST_wrapper( const cudaStream_t streams, + const System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, cuDoubleComplex * d_sol ); -void middle_stuff_DST_wrapper( System sys, +void middle_stuff_DST_wrapper( const cudaStream_t streams, + System sys, const cuDoubleComplex *d_rhat, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y ); -void middle_stuff_ls_DST_wrapper( System sys, +void middle_stuff_ls_DST_wrapper( const cudaStream_t streams, + System sys, const DSTN dst, const cuDoubleComplex *out, const cuDoubleComplex *out2, diff --git a/2D_Program/solver_functions/solver.c b/2D_Program/solver_functions/solver.c index 13f3168..f951123 100644 --- a/2D_Program/solver_functions/solver.c +++ b/2D_Program/solver_functions/solver.c @@ -10,33 +10,14 @@ #include "cuda_kernels.h" void DST( DSTN dst, double _Complex *b, double _Complex *bhat, fftw_plan plan, double *in, fftw_complex *out ); -// void forwardDST( System sys, -// DSTN dst, -// cuDoubleComplex *rhs, -// cuDoubleComplex *bhat, -// fftw_plan plan, -// double * in, -// fftw_complex * out, -// fftw_plan plan2, -// double * in2, -// fftw_complex * out2 ); -// void reverseDST( System sys, -// DSTN dst, -// cuDoubleComplex *xhat, -// cuDoubleComplex *sol, -// fftw_plan plan, -// double * in, -// fftw_complex * out, -// fftw_plan plan2, -// double * in2, -// fftw_complex * out2 ); #ifdef USE_COMBINE void fullDST( const cudaStream_t *streams, + const cudaEvent_t * events, const System sys, const DSTN dst, - const fftw_plan plan, - const fftw_plan plan2, + const cufftHandle plan, + cuDoubleComplex * d_workspace[2], cuDoubleComplex * d_y, double * in, fftw_complex * out, @@ -44,10 +25,10 @@ void fullDST( const cudaStream_t *streams, fftw_complex * out2 ); #else void fullDST( const cudaStream_t *streams, + const cudaEvent_t * events, const System sys, const DSTN dst, - const fftw_plan plan, - const fftw_plan plan2, + const cufftHandle plan, cuDoubleComplex * d_rhat, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y, @@ -62,12 +43,6 @@ void solver( System sys ) { PUSH_RANGE( "solver", 0 ) - int num_streams = 4; - cudaStream_t streams[num_streams]; - for ( int i = 0; i < num_streams; i++ ) { - CUDA_RT_CALL( cudaStreamCreateWithFlags( &streams[i], cudaStreamNonBlocking ) ); - } - DSTN dst; int Nx = sys.lat.Nx, Ny = sys.lat.Ny; //, Nxy = sys.lat.Nxy; @@ -81,20 +56,27 @@ void solver( System sys ) { double * in, *in2; fftw_complex *out, *out2; - CUDA_RT_CALL( cudaMallocManaged( ( void ** )&in, size_in, 1 ) ); - CUDA_RT_CALL( cudaMallocManaged( ( void ** )&in2, size_in, 1 ) ); - CUDA_RT_CALL( cudaMallocManaged( ( void ** )&out, size_out, 1 ) ); - CUDA_RT_CALL( cudaMallocManaged( ( void ** )&out2, size_out, 1 ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( in, size_in, 0, streams[0] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( in2, size_in, 0, streams[0] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( out, size_out, 0, streams[1] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( out2, size_out, 0, streams[1] ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &in ), size_in ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &in2 ), size_in ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &out ), size_out ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &out2 ), size_out ) ); + + PUSH_RANGE( "stream creation", 7 ) + int num_streams = 5; + cudaStream_t streams[num_streams]; + for ( int i = 0; i < num_streams; i++ ) { + CUDA_RT_CALL( cudaStreamCreateWithFlags( &streams[i], cudaStreamNonBlocking ) ); + } + POP_RANGE - CUDA_RT_CALL( cudaMemsetAsync( in, size_in, 0, streams[0] ) ); - CUDA_RT_CALL( cudaMemsetAsync( in2, size_in, 0, streams[0] ) ); - CUDA_RT_CALL( cudaMemsetAsync( out, size_out, 0, streams[1] ) ); - CUDA_RT_CALL( cudaMemsetAsync( out2, size_out, 0, streams[1] ) ); + PUSH_RANGE( "stream creation", 7 ) + int num_events = 5; + cudaEvent_t events[num_events]; + for ( int i = 0; i < num_events; i++ ) { + CUDA_RT_CALL( cudaEventCreateWithFlags( &events[i], cudaEventDisableTiming ) ); + } + POP_RANGE cuDoubleComplex *d_y; CUDA_RT_CALL( cudaMalloc( ( void ** )( &d_y ), sys.lat.Nxy * sizeof( cuDoubleComplex ) ) ); @@ -107,13 +89,19 @@ void solver( System sys ) { CUDA_RT_CALL( cudaMalloc( ( void ** )( &d_xhat ), sys.lat.Nxy * sizeof( cuDoubleComplex ) ) ); #endif - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[0] ) ); + CUDA_RT_CALL( cudaMemsetAsync( in, size_in, 0, NULL ) ); + CUDA_RT_CALL( cudaMemsetAsync( in2, size_in, 0, NULL ) ); + CUDA_RT_CALL( cudaMemsetAsync( out, size_out, 0, NULL ) ); + CUDA_RT_CALL( cudaMemsetAsync( out2, size_out, 0, NULL ) ); + + // The code should be here, BUT there's a bug in cuFFT + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[4] ) ); /**********************BATCHED***************************/ int rank = 1; /* not 2: we are computing 1d transforms */ @@ -127,23 +115,39 @@ void solver( System sys ) { int *onembed = NULL; /**********************BATCHED***************************/ - fftw_plan plan, plan2; /********************* FFTW *********************/ + PUSH_RANGE( "cufft creation", 7 ) +#ifdef USE_COMBINE + cufftHandle plan; + size_t workspace; + cuDoubleComplex *d_workspace[2]; + + CUDA_RT_CALL( cufftCreate( &plan ) ); + CUDA_RT_CALL( cufftSetAutoAllocation( plan, 0 ) ); + CUDA_RT_CALL( cufftMakePlanMany( + plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, howmany, &workspace ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )&d_workspace[0], workspace ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )&d_workspace[1], workspace ) ); +#else + cufftHandle plan; + size_t workspace; + cuDoubleComplex d_workspace; - PUSH_RANGE( "1st fffw_plan", 1 ) - plan = fftw_plan_many_dft_r2c( - rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, FFTW_ESTIMATE ); - plan2 = fftw_plan_many_dft_r2c( - rank, n, howmany, in2, inembed, istride, idist, out2, onembed, ostride, odist, FFTW_ESTIMATE ); + CUDA_RT_CALL( cufftCreate( &plan ) ); + CUDA_RT_CALL( cufftMakePlanMany( + plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, howmany, &workspace ) ); + +#endif POP_RANGE PUSH_RANGE( "DST", 5 ) #ifdef USE_COMBINE - fullDST( streams, sys, dst, plan, plan2, d_y, in, out, in2, out2 ); + fullDST( streams, events, sys, dst, plan, d_workspace, d_y, in, out, in2, out2 ); #else - fullDST( streams, sys, dst, plan, plan2, d_rhat, d_xhat, d_y, in, out, in2, out2 ); + fullDST( streams, events, sys, dst, plan, d_rhat, d_xhat, d_y, in, out, in2, out2 ); #endif - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), cudaCpuDeviceId, NULL ) ); + CUDA_RT_CALL( + cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), cudaCpuDeviceId, streams[0] ) ); POP_RANGE PUSH_RANGE( "Cleanup", 6 ) @@ -162,8 +166,9 @@ void solver( System sys ) { CUDA_RT_CALL( cudaFree( d_xhat ) ); #endif - fftw_destroy_plan( plan ); /********************* FFTW *********************/ - fftw_destroy_plan( plan2 ); /********************* FFTW *********************/ + // for ( int i = 0; i < 2; i++ ) { + CUDA_RT_CALL( cufftDestroy( plan ) ); /********************* FFTW *********************/ + // } POP_RANGE POP_RANGE From 9a88c978cc4ebae5d583bde4afcbf253786fc073 Mon Sep 17 00:00:00 2001 From: Matthew Nicely Date: Tue, 11 May 2021 17:17:25 -0400 Subject: [PATCH 2/3] Move memory push --- 2D_Program/solver_functions/solver.c | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/2D_Program/solver_functions/solver.c b/2D_Program/solver_functions/solver.c index f951123..848cf69 100644 --- a/2D_Program/solver_functions/solver.c +++ b/2D_Program/solver_functions/solver.c @@ -94,15 +94,6 @@ void solver( System sys ) { CUDA_RT_CALL( cudaMemsetAsync( out, size_out, 0, NULL ) ); CUDA_RT_CALL( cudaMemsetAsync( out2, size_out, 0, NULL ) ); - // The code should be here, BUT there's a bug in cuFFT - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); - - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[4] ) ); - /**********************BATCHED***************************/ int rank = 1; /* not 2: we are computing 1d transforms */ int n[] = { N }; @@ -135,10 +126,18 @@ void solver( System sys ) { CUDA_RT_CALL( cufftCreate( &plan ) ); CUDA_RT_CALL( cufftMakePlanMany( plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, howmany, &workspace ) ); - #endif POP_RANGE + // The code should be here, BUT there's a bug in cuFFT + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[4] ) ); + PUSH_RANGE( "DST", 5 ) #ifdef USE_COMBINE fullDST( streams, events, sys, dst, plan, d_workspace, d_y, in, out, in2, out2 ); From 2d563a7399759910cac5771e9f4aef09be30c20d Mon Sep 17 00:00:00 2001 From: Matthew Nicely Date: Wed, 12 May 2021 08:42:34 -0400 Subject: [PATCH 3/3] Optimized out duplicate FFT --- 2D_Program/solver_functions/DST.c | 71 ++++------------- 2D_Program/solver_functions/cuda_kernels.cu | 62 ++++++--------- 2D_Program/solver_functions/cuda_kernels.h | 10 +-- 2D_Program/solver_functions/solver.c | 86 +++++++++------------ 4 files changed, 77 insertions(+), 152 deletions(-) diff --git a/2D_Program/solver_functions/DST.c b/2D_Program/solver_functions/DST.c index c6cb2f2..efa45e7 100644 --- a/2D_Program/solver_functions/DST.c +++ b/2D_Program/solver_functions/DST.c @@ -12,67 +12,28 @@ void fullDST( const cudaStream_t *streams, const System sys, const DSTN dst, const cufftHandle plan, - cuDoubleComplex * d_workspace[2], cuDoubleComplex * d_y, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ) { + fftw_complex * out ) { PUSH_RANGE( "1st DST", 2 ) - CUDA_RT_CALL( cudaEventRecord( events[2], streams[2] ) ); + CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for sys.rhs - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan creation - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 creation - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait for sys.rhs - - load_1st_DST_wrapper( streams[0], sys, dst, sys.rhs, in, in2 ); - - CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); - - CUDA_RT_CALL( cudaStreamWaitEvent( streams[1], events[0], cudaEventWaitDefault ) ); // Wait for load_1st_DST - CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); - CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[0] ) ); + load_1st_DST_wrapper( streams[0], sys, dst, sys.rhs, in ); CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] - - CUDA_RT_CALL( cufftSetStream( plan, streams[1] ) ); - CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[1] ) ); - CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] - - CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); - CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); - CUDA_RT_CALL( cudaEventRecord( events[2], streams[3] ) ); POP_RANGE PUSH_RANGE( "Trig Solver", 3 ) - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan execution - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 execution CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait forsys.U, sys.L, sys.Up - middle_stuff_ls_DST_wrapper( streams[0], sys, dst, out, out2, in, in2, d_y ); - - CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); + middle_stuff_ls_DST_wrapper( streams[0], sys, dst, out, in, d_y ); POP_RANGE PUSH_RANGE( "2nd DST", 4 ) - CUDA_RT_CALL( cudaStreamWaitEvent( streams[1], events[0], cudaEventWaitDefault ) ); // Wait for middle_stuff_ls_DST - CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); - - CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[0] ) ); CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] - CUDA_RT_CALL( cufftSetStream( plan, streams[1] ) ); - CUDA_RT_CALL( cufftSetWorkArea( plan, d_workspace[1] ) ); - CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] - - CUDA_RT_CALL( cudaEventRecord( events[0], streams[0] ) ); - CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); - CUDA_RT_CALL( cudaEventRecord( events[2], streams[4] ) ); - - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[0], cudaEventWaitDefault ) ); // Wait for plan execution - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for plan2 execution - CUDA_RT_CALL( cudaStreamWaitEvent( streams[0], events[2], cudaEventWaitDefault ) ); // Wait for sys.sol - store_2st_DST_wrapper( streams[0], sys, dst, out, out2, sys.sol ); + CUDA_RT_CALL( cudaStreamWaitEvent( streams[3], events[3], cudaEventWaitDefault ) ); // Wait for sys.sol + store_2st_DST_wrapper( streams[0], sys, dst, out, sys.sol ); CUDA_RT_CALL( cudaStreamSynchronize( streams[0] ) ); // Wait for store_2st_DST POP_RANGE @@ -87,17 +48,13 @@ void fullDST( const cudaStream_t *streams, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ) { + fftw_complex * out ) { PUSH_RANGE( "1st DST", 2 ) - load_1st_DST_wrapper( NULL, sys, dst, sys.rhs, in, in2 ); - - CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] - CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] + load_1st_DST_wrapper( NULL, sys, dst, sys.rhs, in ); - store_1st_DST_wrapper( NULL, sys, dst, out, out2, d_rhat ); + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + store_1st_DST_wrapper( NULL, sys, dst, out, d_rhat ); POP_RANGE PUSH_RANGE( "Trig Solver", 3 ) @@ -105,12 +62,10 @@ void fullDST( const cudaStream_t *streams, POP_RANGE PUSH_RANGE( "2nd DST", 4 ) - load_2st_DST_wrapper( NULL, sys, dst, d_xhat, in, in2 ); - - CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] - CUDA_RT_CALL( cufftExecD2Z( plan, in2, out2 ) ); // Running in streams[1] + load_2st_DST_wrapper( NULL, sys, dst, d_xhat, in ); - store_2st_DST_wrapper( NULL, sys, dst, out, out2, sys.sol ); + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + store_2st_DST_wrapper( NULL, sys, dst, out, sys.sol ); POP_RANGE } #endif diff --git a/2D_Program/solver_functions/cuda_kernels.cu b/2D_Program/solver_functions/cuda_kernels.cu index a4327ac..8158f9e 100644 --- a/2D_Program/solver_functions/cuda_kernels.cu +++ b/2D_Program/solver_functions/cuda_kernels.cu @@ -23,8 +23,7 @@ __global__ void __launch_bounds__( 256 ) load_1st_DST( const int N, const int Nx, const int Ny, const cuDoubleComplex *__restrict__ rhs, - double *__restrict__ in, - double *__restrict__ in2 ) { + double *__restrict__ in ) { const int tx { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; const int strideX { static_cast( blockDim.x * gridDim.x ) }; @@ -33,8 +32,8 @@ __global__ void __launch_bounds__( 256 ) load_1st_DST( const int N, for ( int tidY = ty; tidY < Ny; tidY += strideY ) { for ( int tidX = tx; tidX < Nx; tidX += strideX ) { - in[tidY * N + tidX + 1] = rhs[tidX + tidY * Nx].x; - in2[tidY * N + tidX + 1] = rhs[tidX + tidY * Nx].y; + in[tidY * N + tidX + 1] = rhs[tidX + tidY * Nx].x; + in[( N * Ny ) + ( tidY * N + tidX + 1 )] = rhs[tidX + tidY * Nx].y; } } } @@ -54,7 +53,6 @@ __global__ void __launch_bounds__( 256 ) store_1st_DST( const int N, const int NC, const double coef, const cuDoubleComplex *__restrict__ out, - const cuDoubleComplex *__restrict__ out2, cuDoubleComplex *__restrict__ d_rhat ) { const int tx { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; const int strideX { static_cast( blockDim.x * gridDim.x ) }; @@ -65,7 +63,7 @@ __global__ void __launch_bounds__( 256 ) store_1st_DST( const int N, for ( int tidY = ty; tidY < Ny; tidY += strideY ) { for ( int tidX = tx; tidX < Nx; tidX += strideX ) { d_rhat[Nx * tidY + tidX].x = coef * -out[tidY * NC + tidX + 1].y; - d_rhat[Nx * tidY + tidX].y = coef * -out2[tidY * NC + tidX + 1].y; + d_rhat[Nx * tidY + tidX].y = coef * -out[( NC * Ny ) + ( tidY * NC + tidX + 1 )].y; } } } @@ -83,8 +81,7 @@ __global__ void __launch_bounds__( 256 ) load_2st_DST( const int N, const int Nx, const int Ny, const cuDoubleComplex *__restrict__ xhat, - double *__restrict__ in, - double *__restrict__ in2 ) { + double *__restrict__ in ) { const int tx { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; const int strideX { static_cast( blockDim.x * gridDim.x ) }; @@ -93,8 +90,8 @@ __global__ void __launch_bounds__( 256 ) load_2st_DST( const int N, for ( int tidY = ty; tidY < Ny; tidY += strideY ) { for ( int tidX = tx; tidX < Nx; tidX += strideX ) { - in[tidY * N + tidX + 1] = xhat[tidY + tidX * Ny].x; - in2[tidY * N + tidX + 1] = xhat[tidY + tidX * Ny].y; + in[tidY * N + tidX + 1] = xhat[tidY + tidX * Ny].x; + in[( N * Ny ) + ( tidY * N + tidX + 1 )] = xhat[tidY + tidX * Ny].y; } } } @@ -114,7 +111,6 @@ __global__ void __launch_bounds__( 256 ) store_2st_DST( const int N, const int NC, const double coef, const cuDoubleComplex *__restrict__ out, - const cuDoubleComplex *__restrict__ out2, cuDoubleComplex *__restrict__ d_sol ) { const int tx { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; const int strideX { static_cast( blockDim.x * gridDim.x ) }; @@ -125,7 +121,7 @@ __global__ void __launch_bounds__( 256 ) store_2st_DST( const int N, for ( int tidY = ty; tidY < Ny; tidY += strideY ) { for ( int tidX = tx; tidX < Nx; tidX += strideX ) { d_sol[Nx * tidY + tidX].x = coef * -out[tidY * NC + tidX + 1].y; - d_sol[Nx * tidY + tidX].y = coef * -out2[tidY * NC + tidX + 1].y; + d_sol[Nx * tidY + tidX].y = coef * -out[( NC * Ny ) + ( tidY * NC + tidX + 1 )].y; } } } @@ -211,13 +207,11 @@ __global__ void __launch_bounds__( 64 ) middle_stuff_ls_DST( const int N, const int NC, const double coef, const cuDoubleComplex *__restrict__ out, - const cuDoubleComplex *__restrict__ out2, const cuDoubleComplex *__restrict__ d_SysU, const cuDoubleComplex *__restrict__ d_SysL, const cuDoubleComplex *__restrict__ d_SysUp, cuDoubleComplex *__restrict__ d_y, - double *__restrict__ in, - double *__restrict__ in2 ) { + double *__restrict__ in ) { const int tidX { static_cast( blockIdx.x * blockDim.x + threadIdx.x ) }; @@ -226,29 +220,29 @@ __global__ void __launch_bounds__( 64 ) middle_stuff_ls_DST( const int N, if ( tidX < Nx ) { int mx = Ny * tidX; - temp = make_cuDoubleComplex( -out[tidX + 1].y, -out2[tidX + 1].y ); + temp = make_cuDoubleComplex( -out[tidX + 1].y, -out[( NC * Ny ) + tidX + 1].y ); temp = ComplexScale( temp, coef ); d_y[tidX] = temp; #pragma unroll 8 for ( int j = 1; j < Ny; j++ ) { cuDoubleComplex temp2 = cuCmul( d_SysL[mx + j], d_y[( j - 1 ) * Ny + tidX] ); - temp = make_cuDoubleComplex( -out[j * NC + tidX + 1].y, -out2[j * NC + tidX + 1].y ); - temp = ComplexScale( temp, coef ); - d_y[j * Ny + tidX] = cuCsub( temp, temp2 ); + temp = make_cuDoubleComplex( -out[j * NC + tidX + 1].y, -out[( NC * Ny ) + ( j * NC + tidX + 1 )].y ); + temp = ComplexScale( temp, coef ); + d_y[j * Ny + tidX] = cuCsub( temp, temp2 ); } temp = cuCdiv( d_y[( Ny - 1 ) * Ny + tidX], d_SysU[mx + ( Ny - 1 )] ); - in[( Ny - 1 ) * N + tidX + 1] = temp.x; - in2[( Ny - 1 ) * N + tidX + 1] = temp.y; + in[( Ny - 1 ) * N + tidX + 1] = temp.x; + in[( N * Ny ) + ( ( Ny - 1 ) * N + tidX + 1 )] = temp.y; #pragma unroll 4 for ( int j = Ny - 2; j >= 0; j-- ) { cuDoubleComplex temp2 = cuCdiv( cuCsub( d_y[j * Ny + tidX], cuCmul( d_SysUp[mx + j], temp ) ), d_SysU[mx + j] ); - in[j * N + tidX + 1] = temp2.x; - in2[j * N + tidX + 1] = temp2.y; - temp = temp2; + in[j * N + tidX + 1] = temp2.x; + in[( N * Ny ) + ( j * N + tidX + 1 )] = temp2.y; + temp = temp2; } } } @@ -257,8 +251,7 @@ void load_1st_DST_wrapper( const cudaStream_t stream, const System sys, const DSTN dst, const cuDoubleComplex *d_rhs, - double * in, - double * in2 ) { + double * in ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -269,7 +262,7 @@ void load_1st_DST_wrapper( const cudaStream_t stream, dim3 threadPerBlock { 16, 16 }; dim3 blocksPerGrid( numSMs, numSMs ); - void *args[] { &N, &Nx, &Ny, &d_rhs, &in, &in2 }; + void *args[] { &N, &Nx, &Ny, &d_rhs, &in }; CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); @@ -280,7 +273,6 @@ void store_1st_DST_wrapper( const cudaStream_t stream, const System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, cuDoubleComplex * d_rhat ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -294,7 +286,7 @@ void store_1st_DST_wrapper( const cudaStream_t stream, double coef = dst.coef; - void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &d_rhat }; + void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &d_rhat }; CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_1st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); @@ -305,8 +297,7 @@ void load_2st_DST_wrapper( const cudaStream_t stream, const System sys, const DSTN dst, const cuDoubleComplex *d_xhat, - double * in, - double * in2 ) { + double * in ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -317,7 +308,7 @@ void load_2st_DST_wrapper( const cudaStream_t stream, dim3 threadPerBlock { 16, 16 }; dim3 blocksPerGrid( numSMs, numSMs ); - void *args[] { &N, &Nx, &Ny, &d_xhat, &in, &in2 }; + void *args[] { &N, &Nx, &Ny, &d_xhat, &in }; CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &load_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); @@ -328,7 +319,6 @@ void store_2st_DST_wrapper( const cudaStream_t stream, const System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, cuDoubleComplex * d_sol ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -342,7 +332,7 @@ void store_2st_DST_wrapper( const cudaStream_t stream, double coef = dst.coef; - void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &d_sol }; + void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &d_sol }; CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &store_2st_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); @@ -375,9 +365,7 @@ void middle_stuff_ls_DST_wrapper( const cudaStream_t stream, System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, double * in, - double * in2, cuDoubleComplex * d_y ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -391,7 +379,7 @@ void middle_stuff_ls_DST_wrapper( const cudaStream_t stream, int threadPerBlock { 64 }; int blocksPerGrid { ( N + threadPerBlock - 1 ) / threadPerBlock }; - void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &out2, &sys.U, &sys.L, &sys.Up, &d_y, &in, &in2 }; + void *args[] { &N, &Nx, &Ny, &NC, &coef, &out, &sys.U, &sys.L, &sys.Up, &d_y, &in }; CUDA_RT_CALL( cudaLaunchKernel( ( void * )( &middle_stuff_ls_DST ), blocksPerGrid, threadPerBlock, args, 0, stream ) ); diff --git a/2D_Program/solver_functions/cuda_kernels.h b/2D_Program/solver_functions/cuda_kernels.h index 6715ebf..4b93bb0 100644 --- a/2D_Program/solver_functions/cuda_kernels.h +++ b/2D_Program/solver_functions/cuda_kernels.h @@ -10,25 +10,21 @@ void load_1st_DST_wrapper( const cudaStream_t streams, const System sys, const DSTN dst, const cuDoubleComplex *d_rhs, - double * in, - double * in2 ); + double * in ); void store_1st_DST_wrapper( const cudaStream_t streams, const System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, cuDoubleComplex * d_rhat ); void load_2st_DST_wrapper( const cudaStream_t streams, const System sys, const DSTN dst, const cuDoubleComplex *d_xhat, - double * in, - double * in2 ); + double * in ); void store_2st_DST_wrapper( const cudaStream_t streams, const System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, cuDoubleComplex * d_sol ); void middle_stuff_DST_wrapper( const cudaStream_t streams, System sys, @@ -39,9 +35,7 @@ void middle_stuff_ls_DST_wrapper( const cudaStream_t streams, System sys, const DSTN dst, const cuDoubleComplex *out, - const cuDoubleComplex *out2, double * in, - double * in2, cuDoubleComplex * d_y ); #ifdef __cplusplus diff --git a/2D_Program/solver_functions/solver.c b/2D_Program/solver_functions/solver.c index 848cf69..078428c 100644 --- a/2D_Program/solver_functions/solver.c +++ b/2D_Program/solver_functions/solver.c @@ -17,12 +17,9 @@ void fullDST( const cudaStream_t *streams, const System sys, const DSTN dst, const cufftHandle plan, - cuDoubleComplex * d_workspace[2], cuDoubleComplex * d_y, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ); + fftw_complex * out ); #else void fullDST( const cudaStream_t *streams, const cudaEvent_t * events, @@ -33,9 +30,7 @@ void fullDST( const cudaStream_t *streams, cuDoubleComplex * d_xhat, cuDoubleComplex * d_y, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ); + fftw_complex * out ); #endif #if USE_CUFFTW @@ -54,13 +49,11 @@ void solver( System sys ) { size_t size_in = sizeof( double ) * N * Ny; size_t size_out = sizeof( fftw_complex ) * NC * Ny; - double * in, *in2; - fftw_complex *out, *out2; + double * in; + fftw_complex *out; - CUDA_RT_CALL( cudaMalloc( ( void ** )( &in ), size_in ) ); - CUDA_RT_CALL( cudaMalloc( ( void ** )( &in2 ), size_in ) ); - CUDA_RT_CALL( cudaMalloc( ( void ** )( &out ), size_out ) ); - CUDA_RT_CALL( cudaMalloc( ( void ** )( &out2 ), size_out ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &in ), size_in * 2 ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &out ), size_out * 2 ) ); PUSH_RANGE( "stream creation", 7 ) int num_streams = 5; @@ -89,10 +82,8 @@ void solver( System sys ) { CUDA_RT_CALL( cudaMalloc( ( void ** )( &d_xhat ), sys.lat.Nxy * sizeof( cuDoubleComplex ) ) ); #endif - CUDA_RT_CALL( cudaMemsetAsync( in, size_in, 0, NULL ) ); - CUDA_RT_CALL( cudaMemsetAsync( in2, size_in, 0, NULL ) ); - CUDA_RT_CALL( cudaMemsetAsync( out, size_out, 0, NULL ) ); - CUDA_RT_CALL( cudaMemsetAsync( out2, size_out, 0, NULL ) ); + CUDA_RT_CALL( cudaMemsetAsync( in, size_in, 0, streams[0] ) ); + CUDA_RT_CALL( cudaMemsetAsync( out, size_out, 0, streams[1] ) ); /**********************BATCHED***************************/ int rank = 1; /* not 2: we are computing 1d transforms */ @@ -107,54 +98,53 @@ void solver( System sys ) { /**********************BATCHED***************************/ PUSH_RANGE( "cufft creation", 7 ) -#ifdef USE_COMBINE - cufftHandle plan; - size_t workspace; - cuDoubleComplex *d_workspace[2]; + cufftHandle plan; + size_t workspace; CUDA_RT_CALL( cufftCreate( &plan ) ); - CUDA_RT_CALL( cufftSetAutoAllocation( plan, 0 ) ); - CUDA_RT_CALL( cufftMakePlanMany( - plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, howmany, &workspace ) ); - CUDA_RT_CALL( cudaMalloc( ( void ** )&d_workspace[0], workspace ) ); - CUDA_RT_CALL( cudaMalloc( ( void ** )&d_workspace[1], workspace ) ); -#else - cufftHandle plan; - size_t workspace; - cuDoubleComplex d_workspace; - CUDA_RT_CALL( cufftCreate( &plan ) ); +#ifdef USE_COMBINE + CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); CUDA_RT_CALL( cufftMakePlanMany( - plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, howmany, &workspace ) ); -#endif - POP_RANGE + plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, ( howmany * 2 ), &workspace ) ); // The code should be here, BUT there's a bug in cuFFT - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[1] ) ); + CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[2] ) ); + CUDA_RT_CALL( cudaEventRecord( events[2], streams[2] ) ); - CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[4] ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + CUDA_RT_CALL( cudaEventRecord( events[3], streams[3] ) ); PUSH_RANGE( "DST", 5 ) -#ifdef USE_COMBINE - fullDST( streams, events, sys, dst, plan, d_workspace, d_y, in, out, in2, out2 ); -#else - fullDST( streams, events, sys, dst, plan, d_rhat, d_xhat, d_y, in, out, in2, out2 ); -#endif - + fullDST( streams, events, sys, dst, plan, d_y, in, out ); CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), cudaCpuDeviceId, streams[0] ) ); POP_RANGE +#else + CUDA_RT_CALL( cufftMakePlanMany( + plan, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_D2Z, ( howmany * 2 ), &workspace ) ); + + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.rhs, sys.lat.Nxy * sizeof( double _Complex ), 0, NULL ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.U, sys.lat.Nxy * sizeof( double _Complex ), 0, NULL ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.L, sys.lat.Nxy * sizeof( double _Complex ), 0, NULL ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.Up, sys.lat.Nxy * sizeof( double _Complex ), 0, NULL ) ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, NULL ) ); + + PUSH_RANGE( "DST", 5 ) + fullDST( streams, events, sys, dst, plan, d_rhat, d_xhat, d_y, in, out ); + CUDA_RT_CALL( cudaMemPrefetchAsync( sys.sol, sys.lat.Nxy * sizeof( double _Complex ), cudaCpuDeviceId, NULL ) ); + POP_RANGE +#endif PUSH_RANGE( "Cleanup", 6 ) CUDA_RT_CALL( cudaFree( in ) ); CUDA_RT_CALL( cudaFree( out ) ); - CUDA_RT_CALL( cudaFree( in2 ) ); - CUDA_RT_CALL( cudaFree( out2 ) ); for ( int i = 0; i < num_streams; i++ ) { CUDA_RT_CALL( cudaStreamDestroy( streams[i] ) ); @@ -165,9 +155,7 @@ void solver( System sys ) { CUDA_RT_CALL( cudaFree( d_xhat ) ); #endif - // for ( int i = 0; i < 2; i++ ) { CUDA_RT_CALL( cufftDestroy( plan ) ); /********************* FFTW *********************/ - // } POP_RANGE POP_RANGE