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..efa45e7 100644 --- a/2D_Program/solver_functions/DST.c +++ b/2D_Program/solver_functions/DST.c @@ -8,74 +8,64 @@ #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_y, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ) { + fftw_complex * out ) { - 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( cudaStreamWaitEvent( streams[0], events[1], cudaEventWaitDefault ) ); // Wait for sys.rhs - CUDA_RT_CALL( cudaStreamSynchronize( streams[1] ) ); - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ + load_1st_DST_wrapper( streams[0], sys, dst, sys.rhs, in ); + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] 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[2], cudaEventWaitDefault ) ); // Wait forsys.U, sys.L, sys.Up + + middle_stuff_ls_DST_wrapper( streams[0], sys, dst, out, in, d_y ); POP_RANGE - PUSH_RANGE( "forwardDST", 4 ) - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ + PUSH_RANGE( "2nd DST", 4 ) + CUDA_RT_CALL( cufftExecD2Z( plan, in, out ) ); // Running in streams[0] + + 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[3] ) ); - store_2st_DST_wrapper( 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, 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 ); - - CUDA_RT_CALL( cudaStreamSynchronize( streams[1] ) ); - fftw_execute( plan ); /********************* FFTW *********************/ - fftw_execute( plan2 ); /********************* FFTW *********************/ - store_1st_DST_wrapper( sys, dst, out, out2, d_rhat ); + fftw_complex * out ) { + + PUSH_RANGE( "1st DST", 2 ) + load_1st_DST_wrapper( NULL, sys, dst, sys.rhs, in ); + + 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( "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 ); - store_2st_DST_wrapper( 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 eb0d0c3..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; } } } @@ -206,54 +202,56 @@ __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__ d_SysU, + const cuDoubleComplex *__restrict__ d_SysL, + const cuDoubleComplex *__restrict__ d_SysUp, + cuDoubleComplex *__restrict__ d_y, + double *__restrict__ in ) { 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 ); + 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; -#pragma unroll 8 + 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; } } } -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 ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -264,18 +262,17 @@ void load_1st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleCompl 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, 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, cuDoubleComplex * d_rhat ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -289,15 +286,18 @@ void store_1st_DST_wrapper( const System sys, 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, 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 ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; int N = 2 * Nx + 2; @@ -308,18 +308,17 @@ void load_2st_DST_wrapper( const System sys, const DSTN dst, const cuDoubleCompl 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, 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, cuDoubleComplex * d_sol ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -333,15 +332,15 @@ void store_2st_DST_wrapper( const System sys, 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, 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,18 +356,16 @@ 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, double * in, - double * in2, cuDoubleComplex * d_y ) { int Nx = sys.lat.Nx, Ny = sys.lat.Ny; @@ -380,13 +377,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 }; + 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, 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..4b93bb0 100644 --- a/2D_Program/solver_functions/cuda_kernels.h +++ b/2D_Program/solver_functions/cuda_kernels.h @@ -6,28 +6,36 @@ 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 ); +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 ); +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, 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 13f3168..078428c 100644 --- a/2D_Program/solver_functions/solver.c +++ b/2D_Program/solver_functions/solver.c @@ -10,51 +10,27 @@ #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_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, 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, double * in, - fftw_complex * out, - double * in2, - fftw_complex * out2 ); + fftw_complex * out ); #endif #if USE_CUFFTW @@ -62,12 +38,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; @@ -79,22 +49,27 @@ 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; - 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 ) ); + double * in; + fftw_complex *out; - 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 * 2 ) ); + CUDA_RT_CALL( cudaMalloc( ( void ** )( &out ), size_out * 2 ) ); - 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_streams = 5; + cudaStream_t streams[num_streams]; + for ( int i = 0; i < num_streams; i++ ) { + CUDA_RT_CALL( cudaStreamCreateWithFlags( &streams[i], cudaStreamNonBlocking ) ); + } + POP_RANGE + + 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 +82,8 @@ 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( 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.sol, sys.lat.Nxy * sizeof( double _Complex ), 0, streams[3] ) ); + 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 */ @@ -127,31 +97,54 @@ void solver( System sys ) { int *onembed = NULL; /**********************BATCHED***************************/ - fftw_plan plan, plan2; /********************* FFTW *********************/ + PUSH_RANGE( "cufft creation", 7 ) + cufftHandle plan; + size_t 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 ); - POP_RANGE + CUDA_RT_CALL( cufftCreate( &plan ) ); - PUSH_RANGE( "DST", 5 ) #ifdef USE_COMBINE - fullDST( streams, sys, dst, plan, plan2, d_y, in, out, in2, out2 ); + CUDA_RT_CALL( cufftSetStream( plan, streams[0] ) ); + CUDA_RT_CALL( cufftMakePlanMany( + 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[1] ) ); + CUDA_RT_CALL( cudaEventRecord( events[1], streams[1] ) ); + + 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[3] ) ); + CUDA_RT_CALL( cudaEventRecord( events[3], streams[3] ) ); + + PUSH_RANGE( "DST", 5 ) + 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 - fullDST( streams, sys, dst, plan, plan2, d_rhat, d_xhat, d_y, in, out, in2, out2 ); -#endif + 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] ) ); @@ -162,8 +155,7 @@ void solver( System sys ) { CUDA_RT_CALL( cudaFree( d_xhat ) ); #endif - fftw_destroy_plan( plan ); /********************* FFTW *********************/ - fftw_destroy_plan( plan2 ); /********************* FFTW *********************/ + CUDA_RT_CALL( cufftDestroy( plan ) ); /********************* FFTW *********************/ POP_RANGE POP_RANGE