yum-mirror/VkFFT

git clone https://git.yummers.dev/yum-mirror/VkFFT

Tolmachev DmitriiMetal support in VkFFTba7001c

develop
7.6 KiB211 linesraw
1//general parts
2#include <stdio.h>
3#include <vector>
4#include <memory>
5#include <string.h>
6#include <chrono>
7#include <thread>
8#include <iostream>
9#ifndef __STDC_FORMAT_MACROS
10#define __STDC_FORMAT_MACROS
11#endif
12#include <inttypes.h>
13
14//CUDA parts
15#include "cuda_runtime.h"
16#include "device_launch_parameters.h"
17#include <cufft.h>
18
19#include "user_benchmark_cuFFT.h"
20
21void user_benchmark_cuFFT(bool file_output, FILE* output, cuFFTUserSystemParameters* userParams, int device_id)
22{
23	
24	cudaSetDevice(device_id);
25	const int num_runs = 3;
26	double benchmark_result[2] = { 0,0 };//averaged result = sum(system_size/iteration_time)/num_benchmark_samples
27	uint64_t storageComplexSize=8;
28	switch (userParams->P) {
29	case 0:
30		storageComplexSize = (2 * sizeof(float));
31		break;
32	case 1:
33		storageComplexSize = (2 * sizeof(double));
34		break;
35	case 2:
36		storageComplexSize = (2 * 2);
37		break;
38	default:
39		storageComplexSize = (2 * sizeof(float));
40		break;
41	}
42    uint64_t bufferSize = 0;
43    if (userParams->R2C) {
44        bufferSize = (uint64_t)(storageComplexSize / 2) * (userParams->X + 2) * userParams->Y * userParams->Z * userParams->B;
45    }
46    else {
47		bufferSize = (uint64_t)storageComplexSize * userParams->X * userParams->Y * userParams->Z * userParams->B;
48    }
49    
50    float* buffer_input = (float*)malloc(bufferSize);
51    for (uint64_t i = 0; i < bufferSize/sizeof(float); i++) {
52        buffer_input[i] = (float)(2 * ((float)rand()) / RAND_MAX - 1.0);
53    }
54	for (int n = 0; n < 2; n++) {
55		double run_time[num_runs][2];
56		for (int r = 0; r < num_runs; r++) {
57			cufftHandle plan;
58			cufftHandle plan2;
59			void* dataC;
60
61			int dims[3];
62			int FFTdim = 1;
63			if (userParams->Y > 1) FFTdim++;
64			if (userParams->Z > 1) FFTdim++;
65			switch (FFTdim) {
66			case 1:
67				dims[0] = userParams->X;
68				dims[1] = 1;
69				dims[2] = 1;
70				break;
71			case 2:
72				dims[0] = userParams->Y;
73				dims[1] = userParams->X;
74				dims[2] = 1;
75				break;
76			case 3:
77				dims[0] = userParams->Z;
78				dims[1] = userParams->Y;
79				dims[2] = userParams->X;
80				break;
81			}
82			
83			cudaMalloc((void**)&dataC, bufferSize);
84			if (cudaGetLastError() != cudaSuccess) {
85				fprintf(stderr, "Cuda error: Failed to allocate\n");
86				return;
87			}
88			cudaMemcpy(dataC, buffer_input, bufferSize, cudaMemcpyHostToDevice);
89			
90			//forward + inverse
91			int iembed[2][3];
92			int istride[2] = { 1, 1 };
93			int idist[2] = {bufferSize / userParams->B / storageComplexSize, bufferSize / userParams->B / storageComplexSize};
94			if (userParams->R2C) idist[0] *= 2;
95			int oembed[2][3];
96			int ostride[2] = { 1, 1 };
97			int odist[2] = { bufferSize / userParams->B / storageComplexSize, bufferSize / userParams->B / storageComplexSize };
98			if (userParams->R2C) odist[1] *= 2;
99			switch (FFTdim) {
100			case 1:
101				iembed[0][0] = (userParams->R2C) ? dims[0] + 2 : dims[0];
102				oembed[0][0] = (userParams->R2C) ? (dims[0] + 2) / 2 : dims[0];
103
104				iembed[1][0] = (userParams->R2C) ? (dims[0] + 2) / 2 : dims[0];
105				oembed[1][0] = (userParams->R2C) ? dims[0] + 2 : dims[0];
106				break;
107			case 2:
108				iembed[0][0] = dims[0];
109				iembed[0][1] = (userParams->R2C) ? dims[1] + 2 : dims[1];
110				oembed[0][0] = dims[0];
111				oembed[0][1] = (userParams->R2C) ? (dims[1] + 2) / 2 : dims[1];
112
113				iembed[1][0] = dims[0];
114				iembed[1][1] = (userParams->R2C) ? (dims[1] + 2) / 2 : dims[1];
115				oembed[1][0] = dims[0];
116				oembed[1][1] = (userParams->R2C) ? dims[1] + 2 : dims[1];
117				break;
118			case 3:
119				iembed[0][0] = idist[0];
120				iembed[0][1] = dims[1];
121				iembed[0][2] = (userParams->R2C) ? dims[2] + 2 : dims[2];
122				oembed[0][0] = odist[0];
123				oembed[0][1] = dims[1];
124				oembed[0][2] = (userParams->R2C) ? (dims[2] + 2)/2 : dims[2];
125
126				iembed[1][0] = idist[0];
127				iembed[1][1] = dims[1];
128				iembed[1][2] = (userParams->R2C) ? (dims[2] + 2)/2 : dims[2];
129				oembed[1][0] = odist[0];
130				oembed[1][1] = dims[1];
131				oembed[1][2] = (userParams->R2C) ? dims[2] + 2 : dims[2];
132				break;
133			}
134			switch (userParams->P) {
135			case 0:
136				if (userParams->R2C) {
137					cufftPlanMany(&plan, FFTdim, dims, iembed[0], istride[0], idist[0], oembed[0], ostride[0], odist[0], CUFFT_R2C, userParams->B);
138					cufftPlanMany(&plan2, FFTdim, dims, iembed[1], istride[1], idist[1], oembed[1], ostride[1], odist[1], CUFFT_C2R, userParams->B);
139				}
140				else {
141					cufftPlanMany(&plan, FFTdim, dims, iembed[0], istride[0], idist[0], oembed[0], ostride[0], odist[0], CUFFT_C2C, userParams->B);
142				}
143				break;
144			case 1:
145				if (userParams->R2C) {
146					cufftPlanMany(&plan, FFTdim, dims, iembed[0], istride[0], idist[0], oembed[0], ostride[0], odist[0], CUFFT_D2Z, userParams->B);
147					cufftPlanMany(&plan2, FFTdim, dims, iembed[1], istride[1], idist[1], oembed[1], ostride[1], odist[1], CUFFT_Z2D, userParams->B);
148				}
149				else
150					cufftPlanMany(&plan, FFTdim, dims, iembed[0], istride[0], idist[0], oembed[0], ostride[0], odist[0], CUFFT_Z2Z, userParams->B);
151				break;
152			}
153
154			float totTime = 0;
155			std::chrono::steady_clock::time_point timeSubmit = std::chrono::steady_clock::now();
156			for (int i = 0; i < userParams->N; i++) {
157				switch (userParams->P) {
158				case 0:
159					if (userParams->R2C){
160						cufftExecR2C(plan, (cufftReal*) dataC, (cufftComplex*) dataC);
161						cufftExecC2R(plan2, (cufftComplex*) dataC, (cufftReal*) dataC);
162					}else{
163						cufftExecC2C(plan, (cufftComplex*) dataC, (cufftComplex*) dataC, -1);
164						cufftExecC2C(plan, (cufftComplex*) dataC, (cufftComplex*) dataC, 1);
165					}
166					break;
167				case 1:
168					if (userParams->R2C){
169						cufftExecD2Z(plan, (cufftDoubleReal*) dataC, (cufftDoubleComplex*) dataC);
170						cufftExecZ2D(plan2, (cufftDoubleComplex*) dataC, (cufftDoubleReal*) dataC);
171					}else{
172						cufftExecZ2Z(plan, (cufftDoubleComplex*) dataC, (cufftDoubleComplex*) dataC, -1);
173						cufftExecZ2Z(plan, (cufftDoubleComplex*) dataC, (cufftDoubleComplex*) dataC, 1);
174					}
175					break;
176				}
177			}
178			cudaDeviceSynchronize();
179			std::chrono::steady_clock::time_point timeEnd = std::chrono::steady_clock::now();
180			totTime = (std::chrono::duration_cast<std::chrono::microseconds>(timeEnd - timeSubmit).count() * 0.001) / userParams->N;
181				
182			run_time[r][0] = totTime;
183			if (n > 0) {
184				if (r == num_runs - 1) {
185					double std_error = 0;
186					double avg_time = 0;
187					for (uint64_t t = 0; t < num_runs; t++) {
188						avg_time += run_time[t][0];
189					}
190					avg_time /= num_runs;
191					for (uint64_t t = 0; t < num_runs; t++) {
192						std_error += (run_time[t][0] - avg_time) * (run_time[t][0] - avg_time);
193					}
194					std_error = sqrt(std_error / num_runs);
195					if (file_output)
196						fprintf(output, "cuFFT System: %" PRIu64 "x%" PRIu64 "x%" PRIu64 " Batch: %" PRIu64 " Buffer: %" PRIu64 " MB avg_time_per_step: %0.3f ms std_error: %0.3f num_iter: %" PRIu64 " benchmark: %" PRIu64 " scaled bandwidth: %0.1f\n",  userParams->X, userParams->Y, userParams->Z, userParams->B, bufferSize / 1024 / 1024, avg_time, std_error, userParams->N, (uint64_t)(((double)bufferSize / 1024) / avg_time), ((double)bufferSize / 1024.0 / 1024.0 / 1.024 * 4 * FFTdim / avg_time));
197
198					printf("cuFFT System: %" PRIu64 "x%" PRIu64 "x%" PRIu64 " Batch: %" PRIu64 " Buffer: %" PRIu64 " MB avg_time_per_step: %0.3f ms std_error: %0.3f num_iter: %" PRIu64 " benchmark: %" PRIu64 " scaled bandwidth: %0.1f\n", userParams->X, userParams->Y, userParams->Z, userParams->B, bufferSize / 1024 / 1024, avg_time, std_error, userParams->N, (uint64_t)(((double)bufferSize / 1024) / avg_time), ((double)bufferSize / 1024.0 / 1024.0 / 1.024 * 4 * FFTdim / avg_time));
199					benchmark_result[0] += ((double)bufferSize / 1024) / avg_time;
200				}
201
202			}
203			cufftDestroy(plan);
204			if (userParams->R2C)
205				cufftDestroy(plan2);
206			cudaFree(dataC);
207			cudaDeviceSynchronize();
208		}
209	}
210	free(buffer_input);
211}