在 GPU 线程内计算小型复数矩阵的行列式。

huangapple
go评论93阅读模式
英文:

Determinant of a small complex matrix within a GPU thread

问题

问题是计算复杂矩阵（5x5、4x4和3x3）的行列式，要求在特定线程内的设备函数中完成。效率至关重要，因为原问题涉及在每个线程的循环内计算数千个这样的矩阵。

这里是使用直接展开（参数为0运行）和余子式展开（参数为1运行）方法的工作示例：

  1. #include <stdio.h>
  2. #include <complex.h>
  3. #include <cstring>
  5. #include <cuda.h>
  6. #include <cuComplex.h>
  7. #include <cuda_runtime.h>
  8. #include <device_launch_parameters.h>
  10. #define N 5
  12. // ...（略过一部分代码）
  14. __host__ __device__ static __inline__ cuDoubleComplex cuCfms(cuDoubleComplex x,
  15.                                                              cuDoubleComplex y,
  16.                                                              cuDoubleComplex d){
  17.     // ...（略过一部分代码）
  18. }
  20. __host__ __device__ static __inline__ cuDoubleComplex cuCaaaa(cuDoubleComplex x,
  21.                                                               cuDoubleComplex y,
  22.                                                               cuDoubleComplex u,
  23.                                                               cuDoubleComplex v) {
  24.     // ...（略过一部分代码）
  25. }
  27. // ...（略过一部分代码）
  29. __device__ cuDoubleComplex direct_det5 (cuDoubleComplex *B){
  30.     // ...（略过一部分代码）
  31. }
  33. __device__ cuDoubleComplex cuC_5pos (cuDoubleComplex a, cuDoubleComplex b,
  34.                                    cuDoubleComplex c, cuDoubleComplex d,
  35.                                    cuDoubleComplex e) {
  36.     // ...（略过一部分代码）
  37. }
  39. __device__ cuDoubleComplex cuC_5neg (cuDoubleComplex a, cuDoubleComplex b,
  40.                                    cuDoubleComplex c, cuDoubleComplex d,
  41.                                    cuDoubleComplex e) {
  42.     // ...（略过一部分代码）
  43. }
  45. // ...（略过一部分代码）
  47. __device__ cuDoubleComplex det_5x5(cuDoubleComplex *matrix) {
  48.     // ...（略过一部分代码）
  49. }
  51. // ...（略过一部分代码）
  53. __device__ cuDoubleComplex gpu_device_kernel(int i, int method) {
  54.     // ...（略过一部分代码）
  55. }
  57. // ...（略过一部分代码）
  59. __global__ void gpu_global_kernel(cuDoubleComplex *d_res, int method) {
  60.     // ...（略过一部分代码）
  61. }
  63. int main(int argc, char *argv[]) {
  64.     // ...（略过一部分代码）
  65. }

然而，我在考虑这些实现是否真的高效。 我正在寻找一种替代方法，可以在设备函数中计算一般矩阵的行列式。 我提供的示例是最简单直接的方法之一，但我正在寻找另一种方法，比如LU/Cholesky分解，可以从设备函数中调用它。 不幸的是，诸如LAPACK之类的包对于CUDA不可用，但也许来自其他一些软件包的类似兼容函数可以在CUDA设备环境中使用。 对于可以在CUDA设备上使用的类似函数的任何建议或示例将不胜感激！

英文:

The problem is to compute the determinant of a complex matrix (5x5, 4x4, and 3x3) from a device function within a particular thread. Efficiency is crucial, as the original problem involves computing thousands of such matrices in a for loop within each thread.

Here is a working example using direct expansion (run with arg: 0) and cofactor expansion (run with arg: 1) methods:

  1. #include <stdio.h>
  2. #include <complex.h>
  3. #include <cstring>
  4. #include <cuda.h>
  5. #include <cuComplex.h>
  6. #include <cuda_runtime.h>
  7. #include <device_launch_parameters.h>
  8. #define N 5
  9. __host__ __device__ static __inline__  cuDoubleComplex cuCfms( cuDoubleComplex x, 
  10. cuDoubleComplex y, 
  11. cuDoubleComplex d){
  12. double real_res;
  13. double imag_res;
  14. real_res = (cuCreal(x) *  cuCreal(y)) - cuCreal(d);
  15. imag_res = (cuCreal(x) *  cuCimag(y)) - cuCimag(d);
  16. real_res = -(cuCimag(x) * cuCimag(y))  + real_res;  
  17. imag_res =  (cuCimag(x) *  cuCreal(y)) + imag_res;     
  18. return make_cuDoubleComplex(real_res, imag_res);
  19. }
  20. __host__ __device__ static __inline__ cuDoubleComplex cuCaaaa(cuDoubleComplex x,
  21. cuDoubleComplex y,
  22. cuDoubleComplex u,
  23. cuDoubleComplex v) {
  24. return make_cuDoubleComplex (cuCreal(x) + cuCreal(y) + cuCreal(u) + cuCreal(v), 
  25. cuCimag(x) + cuCimag(y) + cuCimag(u) + cuCimag(v));
  26. }
  27. __host__ __device__ static __inline__ cuDoubleComplex cuCasas(cuDoubleComplex x,
  28. cuDoubleComplex y,
  29. cuDoubleComplex u,
  30. cuDoubleComplex v) {
  31. return make_cuDoubleComplex (cuCreal(x) - cuCreal(y) + cuCreal(u) - cuCreal(v), 
  32. cuCimag(x) - cuCimag(y) + cuCimag(u) - cuCimag(v));
  33. }
  34. __host__ __device__ static __inline__ cuDoubleComplex cuCasasa(cuDoubleComplex x,
  35. cuDoubleComplex y,
  36. cuDoubleComplex u,
  37. cuDoubleComplex v,
  38. cuDoubleComplex w) {
  39. return make_cuDoubleComplex (cuCreal(x) - cuCreal(y) + cuCreal(u) - cuCreal(v) + cuCreal(w), 
  40. cuCimag(x) - cuCimag(y) + cuCimag(u) - cuCimag(v) + cuCimag(w) );
  41. }
  42. __device__ cuDoubleComplex cuC_5pos (cuDoubleComplex a, cuDoubleComplex b,
  43. cuDoubleComplex c, cuDoubleComplex d,
  44. cuDoubleComplex e) {
  45. // B[12]*(B[18]*B[24] - B[19]*B[23])
  46. //a*b*c -d*e
  47. return cuCmul(a, cuCfms(b,c, cuCmul(d,e)));
  48. }
  49. __device__ cuDoubleComplex cuC_5neg (cuDoubleComplex a, cuDoubleComplex b,
  50. cuDoubleComplex c, cuDoubleComplex d,
  51. cuDoubleComplex e) {
  52. // B[12]*(B[18]*B[24] - B[19]*B[23])
  53. //a*b*c -d*e
  54. return cuCmul(cuCmul(make_cuDoubleComplex(-1.0,0.0),a), cuCfms(b,c, cuCmul(d,e)));
  55. }
  56. __device__ cuDoubleComplex direct_det5 (cuDoubleComplex *B){
  57. cuDoubleComplex det;
  58. cuDoubleComplex det1 = make_cuDoubleComplex(0.0,0.0);
  59. cuDoubleComplex det2 = make_cuDoubleComplex(0.0,0.0);
  60. cuDoubleComplex det3 = make_cuDoubleComplex(0.0,0.0);
  61. cuDoubleComplex det4 = make_cuDoubleComplex(0.0,0.0);
  62. cuDoubleComplex det5 = make_cuDoubleComplex(0.0,0.0);
  63. cuDoubleComplex det11 = cuCmul(B[0],cuCmul(B[6],cuCadd( cuC_5pos(B[12],B[18],B[24],B[19],B[23]), 
  64. cuCadd(  cuC_5neg(B[13],B[17],B[24],B[22],B[19]), 
  65. cuC_5pos(B[14],B[17],B[23],B[22],B[18]) ) ) ) );
  66. cuDoubleComplex det12 = cuCmul(B[0],cuCmul(B[11],cuCadd(cuC_5pos(B[7],B[18],B[24],B[19],B[23]), 
  67. cuCadd( cuC_5neg(B[8],B[17],B[24],B[22],B[19]), 
  68. cuC_5pos(B[9],B[17],B[23],B[22],B[18]) ) ) ) );
  69. cuDoubleComplex det13 = cuCmul(B[0],cuCmul(B[16],cuCadd(cuC_5pos(B[7],B[13],B[24],B[23],B[14]), 
  70. cuCadd( cuC_5neg(B[8],B[12],B[24],B[22],B[14]), 
  71. cuC_5pos(B[9],B[12],B[23],B[22],B[13]) ) ) ) );
  72. cuDoubleComplex det14 = cuCmul(B[0],cuCmul(B[21],cuCadd(cuC_5pos(B[7],B[13],B[19],B[18],B[14]), 
  73. cuCadd( cuC_5neg(B[8],B[12],B[19],B[17],B[14]), 
  74. cuC_5pos(B[9],B[12],B[18],B[13],B[17]) ) ) ) );
  75. det1 = cuCasas(det11, det12, det13, det14);
  76. cuDoubleComplex det21 = cuCmul(B[1*5+0],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[12],B[18],B[24],B[19],B[23]), cuCadd( cuC_5neg(B[13],B[17],B[24],B[22],B[19]), cuC_5pos(B[14],B[17],B[23],B[22],B[18]) ) ) ) );
  77. cuDoubleComplex det22 = cuCmul(B[1*5+0],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[2],B[18],B[24],B[19],B[23]), cuCadd( cuC_5neg(B[3],B[17],B[24],B[22],B[19]), cuC_5pos(B[4],B[17],B[23],B[22],B[18]) ) ) ) );
  78. cuDoubleComplex det23 = cuCmul(B[1*5+0],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[2],B[13],B[24],B[23],B[14]), cuCadd( cuC_5neg(B[3],B[12],B[24],B[22],B[14]), cuC_5pos(B[4],B[12],B[23],B[22],B[13]) ) ) ) );
  79. cuDoubleComplex det24 = cuCmul(B[1*5+0],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[2],B[13],B[19],B[18],B[14]), cuCadd( cuC_5neg(B[3],B[12],B[19],B[17],B[14]), cuC_5pos(B[4],B[12],B[18],B[13],B[17]) ) ) ) );
  80. det2 = cuCasas(det21, det22, det23, det24);
  81. cuDoubleComplex det31 = cuCmul(B[10],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[3*5+3],B[4*5+4],B[3*5+4],B[4*5+3]), 
  82. cuCadd(  cuC_5neg(B[1*5+3],B[3*5+2],B[4*5+4],B[4*5+2],B[3*5+4]), 
  83. cuC_5pos(B[1*5+4],B[3*5+2],B[4*5+3],B[4*5+2],B[3*5+3]) ) ) ) );
  84. cuDoubleComplex det32 = cuCmul(B[10],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[3*5+3],B[4*5+4],B[3*5+4],B[4*5+3]), 
  85. cuCadd( cuC_5neg(B[0*5+3],B[3*5+2],B[4*5+4],B[4*5+2],B[3*5+4]), 
  86. cuC_5pos(B[0*5+4],B[3*5+2],B[4*5+3],B[4*5+2],B[3*5+3]) ) ) ) );
  87. cuDoubleComplex det33 = cuCmul(B[10],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[4*5+4],B[4*5+3],B[1*5+4]), 
  88. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[4*5+4],B[4*5+2],B[1*5+4]), 
  89. cuC_5pos(B[0*5+4],B[1*5+2],B[4*5+3],B[4*5+2],B[1*5+3]) ) ) ) );
  90. cuDoubleComplex det34 = cuCmul(B[10],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[3*5+4],B[3*5+3],B[1*5+4]), 
  91. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[3*5+4],B[3*5+2],B[1*5+4]), 
  92. cuC_5pos(B[0*5+4],B[1*5+2],B[3*5+3],B[1*5+3],B[3*5+2]) ) ) ) );
  93. det3 = cuCasas(det31, det32, det33, det34);
  94. cuDoubleComplex det41 = cuCmul(B[15],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[2*5+3],B[4*5+4],B[2*5+4],B[4*5+3]), 
  95. cuCadd(  cuC_5neg(B[1*5+3],B[2*5+2],B[4*5+4],B[4*5+2],B[2*5+4]), 
  96. cuC_5pos(B[1*5+4],B[2*5+2],B[4*5+3],B[4*5+2],B[2*5+3]) ) ) ) );
  97. cuDoubleComplex det42 = cuCmul(B[15],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[2*5+3],B[4*5+4],B[2*5+4],B[4*5+3]), 
  98. cuCadd( cuC_5neg(B[0*5+3],B[2*5+2],B[4*5+4],B[4*5+2],B[2*5+4]), 
  99. cuC_5pos(B[0*5+4],B[2*5+2],B[4*5+3],B[4*5+2],B[2*5+3]) ) ) ) );
  100. cuDoubleComplex det43 = cuCmul(B[15],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[4*5+4],B[4*5+3],B[1*5+4]), 
  101. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[4*5+4],B[4*5+2],B[1*5+4]), 
  102. cuC_5pos(B[0*5+4],B[1*5+2],B[4*5+3],B[4*5+2],B[1*5+3]) ) ) ) );
  103. cuDoubleComplex det44 = cuCmul(B[15],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[2*5+4],B[2*5+3],B[1*5+4]), 
  104. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[2*5+4],B[2*5+2],B[1*5+4]), 
  105. cuC_5pos(B[0*5+4],B[1*5+2],B[2*5+3],B[1*5+3],B[2*5+2]) ) ) ) );
  106. det4 = cuCasas(det41, det42, det43, det44);
  107. cuDoubleComplex det51 = cuCmul(B[20],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[2*5+3],B[3*5+4],B[2*5+4],B[3*5+3]), 
  108. cuCadd(  cuC_5neg(B[1*5+3],B[2*5+2],B[3*5+4],B[3*5+2],B[2*5+4]), 
  109. cuC_5pos(B[1*5+4],B[2*5+2],B[3*5+3],B[3*5+2],B[2*5+3]) ) ) ) );
  110. cuDoubleComplex det52 = cuCmul(B[20],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[2*5+3],B[3*5+4],B[2*5+4],B[3*5+3]), 
  111. cuCadd( cuC_5neg(B[0*5+3],B[2*5+2],B[3*5+4],B[3*5+2],B[2*5+4]), 
  112. cuC_5pos(B[0*5+4],B[2*5+2],B[3*5+3],B[3*5+2],B[2*5+3]) ) ) ) );
  113. cuDoubleComplex det53 = cuCmul(B[20],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[3*5+4],B[3*5+3],B[1*5+4]), 
  114. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[3*5+4],B[3*5+2],B[1*5+4]), 
  115. cuC_5pos(B[0*5+4],B[1*5+2],B[3*5+3],B[3*5+2],B[1*5+3]) ) ) ) );
  116. cuDoubleComplex det54 = cuCmul(B[20],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[2*5+4],B[2*5+3],B[1*5+4]), 
  117. cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[2*5+4],B[2*5+2],B[1*5+4]), 
  118. cuC_5pos(B[0*5+4],B[1*5+2],B[2*5+3],B[1*5+3],B[2*5+2]) ) ) ) );
  119. det5 = cuCasas(det51, det52, det53, det54);
  120. det = cuCasasa(det1,det2,det3,det4,det5);
  121. return det;
  122. }
  123. #define N 5
  124. __device__ cuDoubleComplex det_4x4(cuDoubleComplex *matrix) {
  125. cuDoubleComplex det = make_cuDoubleComplex(0.0,0.0);
  126. cuDoubleComplex det1 = make_cuDoubleComplex(0.0,0.0);
  127. cuDoubleComplex det2 = make_cuDoubleComplex(0.0,0.0);
  128. cuDoubleComplex det3 = make_cuDoubleComplex(0.0,0.0);
  129. cuDoubleComplex det4 = make_cuDoubleComplex(0.0,0.0);
  130. // Calculate the determinant using cofactor expansions
  131. det1 = cuCmul(matrix[0], 
  132. cuCadd( cuCsub( cuCmul( matrix[5], 
  133. cuCsub( cuCmul(matrix[10],matrix[15]), 
  134. cuCmul(matrix[11],matrix[14]) ))
  135. ,cuCmul(matrix[6], 
  136. cuCsub( cuCmul(matrix[9],matrix[15]),
  137. cuCmul(matrix[11],matrix[13]) ) ) )
  138. ,cuCmul(matrix[7], 
  139. cuCsub( cuCmul(matrix[9],matrix[14]),
  140. cuCmul(matrix[10],matrix[13]) ) ) ) );
  141. det2 = cuCmul(matrix[1], 
  142. cuCadd( cuCsub( cuCmul( matrix[4], 
  143. cuCsub( cuCmul(matrix[10],matrix[15]), 
  144. cuCmul(matrix[11],matrix[14]) ))
  145. ,cuCmul(matrix[6], 
  146. cuCsub( cuCmul(matrix[8],matrix[15]),
  147. cuCmul(matrix[11],matrix[12]) ) ) )
  148. ,cuCmul(matrix[7], 
  149. cuCsub( cuCmul(matrix[8],matrix[14]),
  150. cuCmul(matrix[10],matrix[12]) ) ) ) );
  151. det3 = cuCmul(matrix[2], 
  152. cuCadd( cuCsub( cuCmul( matrix[4], 
  153. cuCsub( cuCmul(matrix[9],matrix[15]), 
  154. cuCmul(matrix[11],matrix[13]) ))
  155. ,cuCmul(matrix[5], 
  156. cuCsub( cuCmul(matrix[8],matrix[15]),
  157. cuCmul(matrix[11],matrix[12]) ) ) )
  158. ,cuCmul(matrix[7], 
  159. cuCsub( cuCmul(matrix[8],matrix[13]),
  160. cuCmul(matrix[9],matrix[12]) ) ) ) );
  161. det4 = cuCmul(matrix[3], 
  162. cuCadd( cuCsub( cuCmul( matrix[4], 
  163. cuCsub( cuCmul(matrix[9],matrix[14]), 
  164. cuCmul(matrix[10],matrix[13]) ))
  165. ,cuCmul(matrix[5], 
  166. cuCsub( cuCmul(matrix[8],matrix[14]),
  167. cuCmul(matrix[10],matrix[12]) ) ) )
  168. ,cuCmul(matrix[6], 
  169. cuCsub( cuCmul(matrix[8],matrix[13]),
  170. cuCmul(matrix[9],matrix[12]) ) ) ) );
  171. det = cuCsub( cuCadd(det1,det3), cuCadd(det2,det4) );
  172. return det;
  173. }
  174. __device__ cuDoubleComplex det_5x5(cuDoubleComplex *matrix) {
  175. cuDoubleComplex det = make_cuDoubleComplex(0.0,0.0);
  176. // Calculate the determinant using cofactor expansions
  177. for (int col = 0; col < N; col++) {
  178. cuDoubleComplex submatrix[N-1][N-1];
  179. int subrow = 0, subcol = 0;
  180. // Create the submatrix by excluding the current row and column
  181. for (int row = 1; row < N; row++) {
  182. for (int c = 0; c < N; c++) {
  183. if (c == col) {
  184. continue;
  185. }
  186. submatrix[subrow][subcol] = matrix[row * N + c];
  187. subcol++;
  188. }
  189. subrow++;
  190. subcol = 0;
  191. }
  192. // Calculate the cofactor
  193. cuDoubleComplex det_sign = make_cuDoubleComplex((col % 2 == 0 ? 1 : -1),0.0);
  194. cuDoubleComplex cofactor = cuCmul( cuCmul(det_sign,matrix[col]), det_4x4(&submatrix[0][0]));
  195. // Accumulate the cofactors to calculate the determinant
  196. det = cuCadd(det, cofactor);
  197. }
  198. return det;
  199. }
  200. __device__ cuDoubleComplex gpu_device_kernel(int i, int method) {
  201. // some random data, depending on `i`, but for simplification say constant data
  202. cuDoubleComplex matrix[5*5] = {
  203. make_cuDoubleComplex(2.0*(i/69120 * 1000), 2.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(7.0, 0.0),
  204. make_cuDoubleComplex(4.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(8.0, 0.0),
  205. make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(3.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0),
  206. make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(2.0, 0.0),
  207. make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(6.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 2.0)
  208. };
  209. cuDoubleComplex cf_sign[5]={make_cuDoubleComplex(1.0,0.0),
  210. make_cuDoubleComplex(-1.0,0.0),
  211. make_cuDoubleComplex(1.0,0.0),
  212. make_cuDoubleComplex(-1.0,0.0),
  213. make_cuDoubleComplex(1.0,0.0),
  214. };
  215. cuDoubleComplex det;
  216. if (method == 0) {
  217. det = direct_det5(matrix);    
  218. } else {
  219. det = det_5x5(matrix);
  220. }
  221. return det;
  222. }
  223. __global__ void gpu_global_kernel(cuDoubleComplex *d_res, int method) {
  224. // some other computations, irrelevant to current problem
  225. uint64_t total = 69120 * 1000;// 10 * 1000 * 1000; // this `total` is what I am referring to thousands 
  226. cuDoubleComplex det_sum, det;
  227. det_sum = make_cuDoubleComplex(0.0,0.0);
  228. for (int i = 0; i < total; i++) {
  229. det = gpu_device_kernel(i, method);
  230. det_sum = cuCadd(det_sum, det);
  231. }
  232. det_sum = cuCmul(det_sum, make_cuDoubleComplex(1.0/total,0.0) );
  233. d_res[0] = det_sum;
  234. return;
  235. }
  236. int main(int argc, char *argv[]) {
  237. if (argc != 2) {
  238. printf("Usage: %s <num1> <num2>\n", argv[0]);
  239. return 1;
  240. }
  241. int method;
  242. method = atoi(argv[1]);
  243. double complex *res;
  244. res = (double complex*)malloc(sizeof(double complex));
  245. memset(res, 0.0,sizeof(double complex));
  246. cuDoubleComplex *d_res;
  247. cudaMalloc(&d_res, sizeof(cuDoubleComplex) );
  248. int *d_method;
  249. cudaMalloc(&d_method, sizeof(int) );
  250. cudaMemcpy(d_method, &method, sizeof(int), cudaMemcpyHostToDevice);
  251. clock_t direct_start, direct_end;
  252. direct_start = clock();
  253. gpu_global_kernel <<<1,1>>>(d_res, *d_method);
  254. cudaDeviceSynchronize();
  255. direct_end = clock();
  256. cudaMemcpy(res, d_res, sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost );
  257. printf("RESULT: %.16E + %.16E i\n", creal(res[0]), cimag(res[0]));
  258. double direct_time = (double)((double)(direct_end - direct_start) / CLOCKS_PER_SEC);
  259. if (method == 0) {
  260. printf("DIRECT EXEC: %.5f s \n", direct_time);
  261. } else {
  262. printf("COFACTOR EXEC: %.5f s \n", direct_time);        
  263. }    
  264. cudaDeviceReset();
  265. return 0;
  266. }

However, I am wondering if these implementation is even efficient.

I am looking for an alternative method to compute the determinant of a general matrix inside a device function. The example I have provided is one of the most straightforward methods, but an alternative approach, such as LU/Cholesky decomposition, that can be called from the device function is what I am seeking. Unfortunately, packages like LAPACK are not available for CUDA, but maybe a similar compatible function from some other packages could work within the CUDA device environment.

Any suggestions or function examples would be greatly appreciated!

答案1

得分: 1

由于您包含了一个C++头文件(<cstring>)，我假设这是一段C++代码，尽管它看起来非常像C语言。

我直接采用您提供的用于扩展的生成代码，通过简单地定义cuDoubleComplex的*、+和-运算符来实现它。我在下面的代码中将此方法称为“generated”。此外，我对您的代码进行了一些轻微的修改。我添加了一个程序参数来指定GPU线程的数量，以确保GPU的充分利用。我还为每种方法使用了单独的内核。方法0是原始方法，方法1是您的直接方法，方法2是我的“generated”方法。

我使用以下命令编译了该代码：nvcc -std=c++14 -O3 -arch=sm_70 main.cu -o main，并在V100 Pcie上针对32000个GPU线程看到以下时间。

配置 时间
原始 145.01487秒
直接 67.07063秒
生成 24.02771秒

以下是代码：

  1. #include <stdio.h>
  2. #include <complex.h>
  3. #include <cstring>
  5. #include <cuda.h>
  6. #include <cuComplex.h>
  7. #include <cuda_runtime.h>
  8. #include <device_launch_parameters.h>
  10. #define N 5
  12. // ... 其余代码 ...

请注意，我仅返回代码的部分翻译。如果您需要更多信息或有其他问题，请随时提出。

英文:

Since you are including a C++ header (<cstring>) I am assuming its a C++ code, even though it looks very C-like.

I took the general generated code that you provided for expansion and implemented it directly as is by simply defining the operators *,+,- for the cuDoubleComplex .

I call this method "generated" in the code below. Additionally, I made some slight modifications to your code. I added a program parameter to specify the number of GPU threads to have a proper GPU utilization. I also use separate kernels for each method. Method 0 is the original method, method 1 is your direct method, method 2 is my "generated" method.

I compiled the code with nvcc -std=c++14 -O3 -arch=sm_70 main.cu -o main and see the following timings on V100 Pcie for 32000 GPU threads.

Configuration Timing
original 145.01487 s
direct 67.07063 s
generated 24.02771 s

The code:

  1. #include <stdio.h>
  2. #include <complex.h>
  3. #include <cstring>
  5. #include <cuda.h>
  6. #include <cuComplex.h>
  7. #include <cuda_runtime.h>
  8. #include <device_launch_parameters.h>
  10. #define N 5
  13. __host__ __device__
  14. bool operator==(cuDoubleComplex l, cuDoubleComplex r){
  15.     return cuCreal(l) == cuCreal(r) && cuCimag(l) == cuCimag(r);
  16. }
  18. __host__ __device__
  19. cuDoubleComplex operator*(cuDoubleComplex l, cuDoubleComplex r){
  20.     return cuCmul(l,r);
  21. }
  23. __host__ __device__
  24. cuDoubleComplex operator+(cuDoubleComplex l, cuDoubleComplex r){
  25.     return cuCadd(l,r);
  26. }
  28. __host__ __device__
  29. cuDoubleComplex operator-(cuDoubleComplex l, cuDoubleComplex r){
  30.     return cuCsub(l,r);
  31. }
  33. __host__ __device__
  34. cuDoubleComplex operator-(cuDoubleComplex l){
  35.     return make_cuDoubleComplex(-cuCreal(l), -cuCimag(l));
  36. }
  39. template<class Matrix>
  40. __device__ 
  41. cuDoubleComplex generated_det_5x5(Matrix B){
  42.     return 
  43.        B[0*5+0]*(B[1*5+1]*(B[2*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[2*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[2*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))-B[2*5+1]*(B[1*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[1*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[1*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))+B[3*5+1]*(B[1*5+2]*(B[2*5+3]*B[4*5+4]-B[4*5+3]*B[2*5+4])-B[1*5+3]*(B[2*5+2]*B[4*5+4]-B[4*5+2]*B[2*5+4])+B[1*5+4]*(B[2*5+2]*B[4*5+3]-B[4*5+2]*B[2*5+3]))-B[4*5+1]*(B[1*5+2]*(B[2*5+3]*B[3*5+4]-B[3*5+3]*B[2*5+4])-B[1*5+3]*(B[2*5+2]*B[3*5+4]-B[3*5+2]*B[2*5+4])+B[1*5+4]*(B[2*5+2]*B[3*5+3]-B[2*5+3]*B[3*5+2])))
  44.      - B[1*5+0]*(B[0*5+1]*(B[2*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[2*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[2*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))-B[2*5+1]*(B[0*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[0*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[0*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))+B[3*5+1]*(B[0*5+2]*(B[2*5+3]*B[4*5+4]-B[4*5+3]*B[2*5+4])-B[0*5+3]*(B[2*5+2]*B[4*5+4]-B[4*5+2]*B[2*5+4])+B[0*5+4]*(B[2*5+2]*B[4*5+3]-B[4*5+2]*B[2*5+3]))-B[4*5+1]*(B[0*5+2]*(B[2*5+3]*B[3*5+4]-B[3*5+3]*B[2*5+4])-B[0*5+3]*(B[2*5+2]*B[3*5+4]-B[3*5+2]*B[2*5+4])+B[0*5+4]*(B[2*5+2]*B[3*5+3]-B[2*5+3]*B[3*5+2])))
  45.      + B[2*5+0]*(B[0*5+1]*(B[1*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[1*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[1*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))-B[1*5+1]*(B[0*5+2]*(B[3*5+3]*B[4*5+4]-B[3*5+4]*B[4*5+3])-B[0*5+3]*(B[3*5+2]*B[4*5+4]-B[4*5+2]*B[3*5+4])+B[0*5+4]*(B[3*5+2]*B[4*5+3]-B[4*5+2]*B[3*5+3]))+B[3*5+1]*(B[0*5+2]*(B[1*5+3]*B[4*5+4]-B[4*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[4*5+4]-B[4*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[4*5+3]-B[4*5+2]*B[1*5+3]))-B[4*5+1]*(B[0*5+2]*(B[1*5+3]*B[3*5+4]-B[3*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[3*5+4]-B[3*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[3*5+3]-B[1*5+3]*B[3*5+2])))
  46.      - B[3*5+0]*(B[0*5+1]*(B[1*5+2]*(B[2*5+3]*B[4*5+4]-B[2*5+4]*B[4*5+3])-B[1*5+3]*(B[2*5+2]*B[4*5+4]-B[4*5+2]*B[2*5+4])+B[1*5+4]*(B[2*5+2]*B[4*5+3]-B[4*5+2]*B[2*5+3]))-B[1*5+1]*(B[0*5+2]*(B[2*5+3]*B[4*5+4]-B[2*5+4]*B[4*5+3])-B[0*5+3]*(B[2*5+2]*B[4*5+4]-B[4*5+2]*B[2*5+4])+B[0*5+4]*(B[2*5+2]*B[4*5+3]-B[4*5+2]*B[2*5+3]))+B[2*5+1]*(B[0*5+2]*(B[1*5+3]*B[4*5+4]-B[4*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[4*5+4]-B[4*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[4*5+3]-B[4*5+2]*B[1*5+3]))-B[4*5+1]*(B[0*5+2]*(B[1*5+3]*B[2*5+4]-B[2*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[2*5+4]-B[2*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[2*5+3]-B[1*5+3]*B[2*5+2])))
  47.      + B[4*5+0]*(B[0*5+1]*(B[1*5+2]*(B[2*5+3]*B[3*5+4]-B[2*5+4]*B[3*5+3])-B[1*5+3]*(B[2*5+2]*B[3*5+4]-B[3*5+2]*B[2*5+4])+B[1*5+4]*(B[2*5+2]*B[3*5+3]-B[3*5+2]*B[2*5+3]))-B[1*5+1]*(B[0*5+2]*(B[2*5+3]*B[3*5+4]-B[2*5+4]*B[3*5+3])-B[0*5+3]*(B[2*5+2]*B[3*5+4]-B[3*5+2]*B[2*5+4])+B[0*5+4]*(B[2*5+2]*B[3*5+3]-B[3*5+2]*B[2*5+3]))+B[2*5+1]*(B[0*5+2]*(B[1*5+3]*B[3*5+4]-B[3*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[3*5+4]-B[3*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[3*5+3]-B[3*5+2]*B[1*5+3]))-B[3*5+1]*(B[0*5+2]*(B[1*5+3]*B[2*5+4]-B[2*5+3]*B[1*5+4])-B[0*5+3]*(B[1*5+2]*B[2*5+4]-B[2*5+2]*B[1*5+4])+B[0*5+4]*(B[1*5+2]*B[2*5+3]-B[1*5+3]*B[2*5+2])));
  48. }
  57. __host__ __device__ static __inline__  cuDoubleComplex cuCfms( cuDoubleComplex x, 
  58.                                                                cuDoubleComplex y, 
  59.                                                                cuDoubleComplex d){
  60.     double real_res;
  61.     double imag_res;
  62.     
  63.     real_res = (cuCreal(x) *  cuCreal(y)) - cuCreal(d);
  64.     imag_res = (cuCreal(x) *  cuCimag(y)) - cuCimag(d);
  65.             
  66.     real_res = -(cuCimag(x) * cuCimag(y))  + real_res;  
  67.     imag_res =  (cuCimag(x) *  cuCreal(y)) + imag_res;     
  68.      
  69.     return make_cuDoubleComplex(real_res, imag_res);
  70. }
  72. __host__ __device__ static __inline__ cuDoubleComplex cuCaaaa(cuDoubleComplex x,
  73.                                                               cuDoubleComplex y,
  74.                                                               cuDoubleComplex u,
  75.                                                               cuDoubleComplex v) {
  77.     return make_cuDoubleComplex (cuCreal(x) + cuCreal(y) + cuCreal(u) + cuCreal(v), 
  78.                                  cuCimag(x) + cuCimag(y) + cuCimag(u) + cuCimag(v));
  79. }
  81. __host__ __device__ static __inline__ cuDoubleComplex cuCasas(cuDoubleComplex x,
  82.                                                               cuDoubleComplex y,
  83.                                                               cuDoubleComplex u,
  84.                                                               cuDoubleComplex v) {
  86.     return make_cuDoubleComplex (cuCreal(x) - cuCreal(y) + cuCreal(u) - cuCreal(v), 
  87.                                  cuCimag(x) - cuCimag(y) + cuCimag(u) - cuCimag(v));
  88. }
  90. __host__ __device__ static __inline__ cuDoubleComplex cuCasasa(cuDoubleComplex x,
  91.                                                                cuDoubleComplex y,
  92.                                                                cuDoubleComplex u,
  93.                                                                cuDoubleComplex v,
  94.                                                                cuDoubleComplex w) {
  96.     return make_cuDoubleComplex (cuCreal(x) - cuCreal(y) + cuCreal(u) - cuCreal(v) + cuCreal(w), 
  97.                                  cuCimag(x) - cuCimag(y) + cuCimag(u) - cuCimag(v) + cuCimag(w) );
  98. }
  101. __device__ cuDoubleComplex cuC_5pos (cuDoubleComplex a, cuDoubleComplex b,
  102.                                    cuDoubleComplex c, cuDoubleComplex d,
  103.                                    cuDoubleComplex e) {
  104.      // B[12]*(B[18]*B[24] - B[19]*B[23])
  105.     //a*b*c -d*e
  106.     return cuCmul(a, cuCfms(b,c, cuCmul(d,e)));
  107. }
  109. __device__ cuDoubleComplex cuC_5neg (cuDoubleComplex a, cuDoubleComplex b,
  110.                                    cuDoubleComplex c, cuDoubleComplex d,
  111.                                    cuDoubleComplex e) {
  112.      // B[12]*(B[18]*B[24] - B[19]*B[23])
  113.     //a*b*c -d*e
  114.     return cuCmul(cuCmul(make_cuDoubleComplex(-1.0,0.0),a), cuCfms(b,c, cuCmul(d,e)));
  115. }
  117. __device__ cuDoubleComplex direct_det5 (cuDoubleComplex *B){
  118.     cuDoubleComplex det;
  120.     cuDoubleComplex det1 = make_cuDoubleComplex(0.0,0.0);
  121.     cuDoubleComplex det2 = make_cuDoubleComplex(0.0,0.0);
  122.     cuDoubleComplex det3 = make_cuDoubleComplex(0.0,0.0);
  123.     cuDoubleComplex det4 = make_cuDoubleComplex(0.0,0.0);
  124.     cuDoubleComplex det5 = make_cuDoubleComplex(0.0,0.0);
  127.     cuDoubleComplex det11 = cuCmul(B[0],cuCmul(B[6],cuCadd( cuC_5pos(B[12],B[18],B[24],B[19],B[23]), 
  128.                                                    cuCadd(  cuC_5neg(B[13],B[17],B[24],B[22],B[19]), 
  129.                                                             cuC_5pos(B[14],B[17],B[23],B[22],B[18]) ) ) ) );
  130.     cuDoubleComplex det12 = cuCmul(B[0],cuCmul(B[11],cuCadd(cuC_5pos(B[7],B[18],B[24],B[19],B[23]), 
  131.                                                     cuCadd( cuC_5neg(B[8],B[17],B[24],B[22],B[19]), 
  132.                                                             cuC_5pos(B[9],B[17],B[23],B[22],B[18]) ) ) ) );
  133.     cuDoubleComplex det13 = cuCmul(B[0],cuCmul(B[16],cuCadd(cuC_5pos(B[7],B[13],B[24],B[23],B[14]), 
  134.                                                     cuCadd( cuC_5neg(B[8],B[12],B[24],B[22],B[14]), 
  135.                                                             cuC_5pos(B[9],B[12],B[23],B[22],B[13]) ) ) ) );
  136.     cuDoubleComplex det14 = cuCmul(B[0],cuCmul(B[21],cuCadd(cuC_5pos(B[7],B[13],B[19],B[18],B[14]), 
  137.                                                     cuCadd( cuC_5neg(B[8],B[12],B[19],B[17],B[14]), 
  138.                                                             cuC_5pos(B[9],B[12],B[18],B[13],B[17]) ) ) ) );
  140.     det1 = cuCasas(det11, det12, det13, det14);
  143.     cuDoubleComplex det21 = cuCmul(B[1*5+0],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[12],B[18],B[24],B[19],B[23]), cuCadd( cuC_5neg(B[13],B[17],B[24],B[22],B[19]), cuC_5pos(B[14],B[17],B[23],B[22],B[18]) ) ) ) );
  144.     cuDoubleComplex det22 = cuCmul(B[1*5+0],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[2],B[18],B[24],B[19],B[23]), cuCadd( cuC_5neg(B[3],B[17],B[24],B[22],B[19]), cuC_5pos(B[4],B[17],B[23],B[22],B[18]) ) ) ) );
  145.     cuDoubleComplex det23 = cuCmul(B[1*5+0],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[2],B[13],B[24],B[23],B[14]), cuCadd( cuC_5neg(B[3],B[12],B[24],B[22],B[14]), cuC_5pos(B[4],B[12],B[23],B[22],B[13]) ) ) ) );
  146.     cuDoubleComplex det24 = cuCmul(B[1*5+0],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[2],B[13],B[19],B[18],B[14]), cuCadd( cuC_5neg(B[3],B[12],B[19],B[17],B[14]), cuC_5pos(B[4],B[12],B[18],B[13],B[17]) ) ) ) );
  148.     det2 = cuCasas(det21, det22, det23, det24);
  151.     cuDoubleComplex det31 = cuCmul(B[10],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[3*5+3],B[4*5+4],B[3*5+4],B[4*5+3]), 
  152.                                                        cuCadd(  cuC_5neg(B[1*5+3],B[3*5+2],B[4*5+4],B[4*5+2],B[3*5+4]), 
  153.                                                                 cuC_5pos(B[1*5+4],B[3*5+2],B[4*5+3],B[4*5+2],B[3*5+3]) ) ) ) );
  154.     cuDoubleComplex det32 = cuCmul(B[10],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[3*5+3],B[4*5+4],B[3*5+4],B[4*5+3]), 
  155.                                                         cuCadd( cuC_5neg(B[0*5+3],B[3*5+2],B[4*5+4],B[4*5+2],B[3*5+4]), 
  156.                                                                 cuC_5pos(B[0*5+4],B[3*5+2],B[4*5+3],B[4*5+2],B[3*5+3]) ) ) ) );
  157.     cuDoubleComplex det33 = cuCmul(B[10],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[4*5+4],B[4*5+3],B[1*5+4]), 
  158.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[4*5+4],B[4*5+2],B[1*5+4]), 
  159.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[4*5+3],B[4*5+2],B[1*5+3]) ) ) ) );
  160.     cuDoubleComplex det34 = cuCmul(B[10],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[3*5+4],B[3*5+3],B[1*5+4]), 
  161.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[3*5+4],B[3*5+2],B[1*5+4]), 
  162.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[3*5+3],B[1*5+3],B[3*5+2]) ) ) ) );
  165.     det3 = cuCasas(det31, det32, det33, det34);
  168.     cuDoubleComplex det41 = cuCmul(B[15],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[2*5+3],B[4*5+4],B[2*5+4],B[4*5+3]), 
  169.                                                        cuCadd(  cuC_5neg(B[1*5+3],B[2*5+2],B[4*5+4],B[4*5+2],B[2*5+4]), 
  170.                                                                 cuC_5pos(B[1*5+4],B[2*5+2],B[4*5+3],B[4*5+2],B[2*5+3]) ) ) ) );
  171.     cuDoubleComplex det42 = cuCmul(B[15],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[2*5+3],B[4*5+4],B[2*5+4],B[4*5+3]), 
  172.                                                         cuCadd( cuC_5neg(B[0*5+3],B[2*5+2],B[4*5+4],B[4*5+2],B[2*5+4]), 
  173.                                                                 cuC_5pos(B[0*5+4],B[2*5+2],B[4*5+3],B[4*5+2],B[2*5+3]) ) ) ) );
  174.     cuDoubleComplex det43 = cuCmul(B[15],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[4*5+4],B[4*5+3],B[1*5+4]), 
  175.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[4*5+4],B[4*5+2],B[1*5+4]), 
  176.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[4*5+3],B[4*5+2],B[1*5+3]) ) ) ) );
  177.     cuDoubleComplex det44 = cuCmul(B[15],cuCmul(B[4*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[2*5+4],B[2*5+3],B[1*5+4]), 
  178.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[2*5+4],B[2*5+2],B[1*5+4]), 
  179.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[2*5+3],B[1*5+3],B[2*5+2]) ) ) ) );
  182.     det4 = cuCasas(det41, det42, det43, det44);
  185.     cuDoubleComplex det51 = cuCmul(B[20],cuCmul(B[0*5+1],cuCadd(cuC_5pos(B[1*5+2],B[2*5+3],B[3*5+4],B[2*5+4],B[3*5+3]), 
  186.                                                        cuCadd(  cuC_5neg(B[1*5+3],B[2*5+2],B[3*5+4],B[3*5+2],B[2*5+4]), 
  187.                                                                 cuC_5pos(B[1*5+4],B[2*5+2],B[3*5+3],B[3*5+2],B[2*5+3]) ) ) ) );
  188.     cuDoubleComplex det52 = cuCmul(B[20],cuCmul(B[1*5+1],cuCadd(cuC_5pos(B[0*5+2],B[2*5+3],B[3*5+4],B[2*5+4],B[3*5+3]), 
  189.                                                         cuCadd( cuC_5neg(B[0*5+3],B[2*5+2],B[3*5+4],B[3*5+2],B[2*5+4]), 
  190.                                                                 cuC_5pos(B[0*5+4],B[2*5+2],B[3*5+3],B[3*5+2],B[2*5+3]) ) ) ) );
  191.     cuDoubleComplex det53 = cuCmul(B[20],cuCmul(B[2*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[3*5+4],B[3*5+3],B[1*5+4]), 
  192.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[3*5+4],B[3*5+2],B[1*5+4]), 
  193.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[3*5+3],B[3*5+2],B[1*5+3]) ) ) ) );
  194.     cuDoubleComplex det54 = cuCmul(B[20],cuCmul(B[3*5+1],cuCadd(cuC_5pos(B[0*5+2],B[1*5+3],B[2*5+4],B[2*5+3],B[1*5+4]), 
  195.                                                         cuCadd( cuC_5neg(B[0*5+3],B[1*5+2],B[2*5+4],B[2*5+2],B[1*5+4]), 
  196.                                                                 cuC_5pos(B[0*5+4],B[1*5+2],B[2*5+3],B[1*5+3],B[2*5+2]) ) ) ) );
  199.     det5 = cuCasas(det51, det52, det53, det54);
  201.     det = cuCasasa(det1,det2,det3,det4,det5);
  203.     return det;
  204. }
  208. __device__ cuDoubleComplex det_4x4(cuDoubleComplex *matrix) {
  209.     cuDoubleComplex det = make_cuDoubleComplex(0.0,0.0);
  210.     cuDoubleComplex det1 = make_cuDoubleComplex(0.0,0.0);
  211.     cuDoubleComplex det2 = make_cuDoubleComplex(0.0,0.0);
  212.     cuDoubleComplex det3 = make_cuDoubleComplex(0.0,0.0);
  213.     cuDoubleComplex det4 = make_cuDoubleComplex(0.0,0.0);
  215.     // Calculate the determinant using cofactor expansions
  216.     det1 = cuCmul(matrix[0], 
  217.                   cuCadd( cuCsub( cuCmul( matrix[5], 
  218.                                           cuCsub( cuCmul(matrix[10],matrix[15]), 
  219.                                                   cuCmul(matrix[11],matrix[14]) ))
  220.                                  ,cuCmul(matrix[6], 
  221.                                          cuCsub( cuCmul(matrix[9],matrix[15]),
  222.                                                  cuCmul(matrix[11],matrix[13]) ) ) )
  223.                          ,cuCmul(matrix[7], 
  224.                                  cuCsub( cuCmul(matrix[9],matrix[14]),
  225.                                          cuCmul(matrix[10],matrix[13]) ) ) ) );
  229.     det2 = cuCmul(matrix[1], 
  230.                   cuCadd( cuCsub( cuCmul( matrix[4], 
  231.                                           cuCsub( cuCmul(matrix[10],matrix[15]), 
  232.                                                   cuCmul(matrix[11],matrix[14]) ))
  233.                                  ,cuCmul(matrix[6], 
  234.                                          cuCsub( cuCmul(matrix[8],matrix[15]),
  235.                                                  cuCmul(matrix[11],matrix[12]) ) ) )
  236.                          ,cuCmul(matrix[7], 
  237.                                  cuCsub( cuCmul(matrix[8],matrix[14]),
  238.                                          cuCmul(matrix[10],matrix[12]) ) ) ) );
  241.     det3 = cuCmul(matrix[2], 
  242.                   cuCadd( cuCsub( cuCmul( matrix[4], 
  243.                                           cuCsub( cuCmul(matrix[9],matrix[15]), 
  244.                                                   cuCmul(matrix[11],matrix[13]) ))
  245.                                  ,cuCmul(matrix[5], 
  246.                                          cuCsub( cuCmul(matrix[8],matrix[15]),
  247.                                                  cuCmul(matrix[11],matrix[12]) ) ) )
  248.                          ,cuCmul(matrix[7], 
  249.                                  cuCsub( cuCmul(matrix[8],matrix[13]),
  250.                                          cuCmul(matrix[9],matrix[12]) ) ) ) );
  251.     
  252.     
  253.     det4 = cuCmul(matrix[3], 
  254.                   cuCadd( cuCsub( cuCmul( matrix[4], 
  255.                                           cuCsub( cuCmul(matrix[9],matrix[14]), 
  256.                                                   cuCmul(matrix[10],matrix[13]) ))
  257.                                  ,cuCmul(matrix[5], 
  258.                                          cuCsub( cuCmul(matrix[8],matrix[14]),
  259.                                                  cuCmul(matrix[10],matrix[12]) ) ) )
  260.                          ,cuCmul(matrix[6], 
  261.                                  cuCsub( cuCmul(matrix[8],matrix[13]),
  262.                                          cuCmul(matrix[9],matrix[12]) ) ) ) );
  263.     
  264.     det = cuCsub( cuCadd(det1,det3), cuCadd(det2,det4) );
  266.     return det;
  267. }
  270. __device__ cuDoubleComplex det_5x5(cuDoubleComplex *matrix) {
  271.     cuDoubleComplex det = make_cuDoubleComplex(0.0,0.0);
  274.     // Calculate the determinant using cofactor expansions
  275.     #pragma unroll
  276.     for (int col = 0; col < N; col++) {
  277.         cuDoubleComplex submatrix[N-1][N-1];
  278.         int subrow = 0, subcol = 0;
  280.         // Create the submatrix by excluding the current row and column
  281.         #pragma unroll
  282.         for (int row = 1; row < N; row++) {
  283.             #pragma unroll
  284.             for (int c = 0; c < N; c++) {
  285.                 if (c == col) {
  286.                     continue;
  287.                 }
  288.                 submatrix[subrow][subcol] = matrix[row * N + c];
  289.                 subcol++;
  290.             }
  291.             subrow++;
  292.             subcol = 0;
  293.         }
  295.         // Calculate the cofactor
  296.         cuDoubleComplex det_sign = make_cuDoubleComplex((col % 2 == 0 ? 1 : -1),0.0);
  297.         cuDoubleComplex cofactor = cuCmul( cuCmul(det_sign,matrix[col]), det_4x4(&submatrix[0][0]));
  300.         // Accumulate the cofactors to calculate the determinant
  301.         det = cuCadd(det, cofactor);
  302.     }
  305.     return det;
  306. }
  308. __device__ cuDoubleComplex gpu_device_kernel_original(int i) {
  311.     // some random data, depending on `i`, but for simplification say constant data
  312.     cuDoubleComplex matrix[5*5] = {
  313.             make_cuDoubleComplex(2.0*(i/69120 * 1000), 2.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(7.0, 0.0),
  314.             make_cuDoubleComplex(4.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(8.0, 0.0),
  315.             make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(3.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0),
  316.             make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(2.0, 0.0),
  317.             make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(6.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 2.0)
  318.     };
  320.     cuDoubleComplex det = det_5x5(matrix);
  322.     return det;
  323. }
  325. __device__ cuDoubleComplex gpu_device_kernel_direct(int i) {
  328.     // some random data, depending on `i`, but for simplification say constant data
  329.     cuDoubleComplex matrix[5*5] = {
  330.             make_cuDoubleComplex(2.0*(i/69120 * 1000), 2.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(7.0, 0.0),
  331.             make_cuDoubleComplex(4.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(8.0, 0.0),
  332.             make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(3.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0),
  333.             make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(2.0, 0.0),
  334.             make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(6.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 2.0)
  335.     };
  337.     cuDoubleComplex det = direct_det5(matrix);    
  339.     return det;
  340. }
  342. __device__ cuDoubleComplex gpu_device_kernel_generated(int i) {
  345.     // some random data, depending on `i`, but for simplification say constant data
  346.     cuDoubleComplex matrix[5*5] = {
  347.             make_cuDoubleComplex(2.0*(i/69120 * 1000), 2.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(7.0, 0.0),
  348.             make_cuDoubleComplex(4.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(8.0, 0.0),
  349.             make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(3.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(2.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0),
  350.             make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(1.0, 0.0), make_cuDoubleComplex(2.0, 0.0),
  351.             make_cuDoubleComplex(5.0, 0.0), make_cuDoubleComplex(6.0*(i/69120 * 1000), 0.0), make_cuDoubleComplex(4.0, 0.0), make_cuDoubleComplex(3.0, 0.0), make_cuDoubleComplex(2.0, 2.0)
  352.     };
  354.     cuDoubleComplex det = generated_det_5x5(matrix);    
  356.     return det;
  357. }
  360. __global__ void gpu_global_kernel_original(cuDoubleComplex *d_res, int numGpuThreads) {
  361.     const int tid = threadIdx.x + blockIdx.x * blockDim.x;
  362.     if(tid < numGpuThreads){
  364.         // some other computations, irrelevant to current problem
  366.         uint64_t total = 69120 * 1000;// 10 * 1000 * 1000; // this `total` is what I am referring to thousands 
  367.         
  368.         cuDoubleComplex det_sum, det;
  369.         det_sum = make_cuDoubleComplex(0.0,0.0);
  370.         for (int i = 0; i < total; i++) {
  371.             det = gpu_device_kernel_original(i);
  372.             det_sum = cuCadd(det_sum, det);
  373.         }
  375.         det_sum = cuCmul(det_sum, make_cuDoubleComplex(1.0/total,0.0) );
  377.         d_res[tid] = det_sum;
  378.     }
  379. }
  382. __global__ void gpu_global_kernel_direct(cuDoubleComplex *d_res, int numGpuThreads) {
  383.     const int tid = threadIdx.x + blockIdx.x * blockDim.x;
  384.     if(tid < numGpuThreads){
  386.         // some other computations, irrelevant to current problem
  388.         uint64_t total = 69120 * 1000;// 10 * 1000 * 1000; // this `total` is what I am referring to thousands 
  389.         
  390.         cuDoubleComplex det_sum, det;
  391.         det_sum = make_cuDoubleComplex(0.0,0.0);
  392.         for (int i = 0; i < total; i++) {
  393.             det = gpu_device_kernel_direct(i);
  394.             det_sum = cuCadd(det_sum, det);
  395.         }
  397.         det_sum = cuCmul(det_sum, make_cuDoubleComplex(1.0/total,0.0) );
  399.         d_res[tid] = det_sum;
  400.     }
  401. }
  404. __global__ void gpu_global_kernel_generated(cuDoubleComplex *d_res, int numGpuThreads) {
  405.     const int tid = threadIdx.x + blockIdx.x * blockDim.x;
  406.     if(tid < numGpuThreads){
  408.         // some other computations, irrelevant to current problem
  410.         uint64_t total = 69120 * 1000;// 10 * 1000 * 1000; // this `total` is what I am referring to thousands 
  411.         
  412.         cuDoubleComplex det_sum, det;
  413.         det_sum = make_cuDoubleComplex(0.0,0.0);
  414.         for (int i = 0; i < total; i++) {
  415.             det = gpu_device_kernel_generated(i);
  416.             det_sum = cuCadd(det_sum, det);
  417.         }
  419.         det_sum = cuCmul(det_sum, make_cuDoubleComplex(1.0/total,0.0) );
  421.         d_res[tid] = det_sum;
  422.     }
  423. }
  425. int main(int argc, char *argv[]) {
  427.     if (argc < 2) {
  428.         printf("Usage: %s <method> [<numGpuThreads>]\n", argv[0]);
  429.         printf("method: 0 = original, 1 = direct, 2 = generated(no explicit cucomplex)\n");
  430.         return 1;
  431.     }
  433.     int method;
  434.     method = atoi(argv[1]);
  436.     int numGpuThreads = 1;
  437.     if(argc > 2){
  438.         numGpuThreads = atoi(argv[2]);
  439.     }
  440.     
  441.     cuDoubleComplex *res;
  442.     res = (cuDoubleComplex*)malloc(sizeof(cuDoubleComplex) * numGpuThreads);
  443.     memset(res, 0.0,sizeof(cuDoubleComplex) * numGpuThreads);
  446.     cuDoubleComplex *d_res;
  447.     cudaMalloc(&d_res, sizeof(cuDoubleComplex) * numGpuThreads );
  449.     clock_t direct_start, direct_end;
  450.     direct_start = clock();
  451.     dim3 block = 128;
  452.     dim3 grid = (numGpuThreads + block.x -1) / block.x;
  453.     switch(method){
  454.         case 0: printf("original\n"); gpu_global_kernel_original<<<grid,block>>>(d_res, numGpuThreads); break;
  455.         case 1: printf("direct\n"); gpu_global_kernel_direct<<<grid,block>>>(d_res, numGpuThreads); break;
  456.         case 2: printf("generated\n"); gpu_global_kernel_generated<<<grid,block>>>(d_res, numGpuThreads); break;
  457.     }
  459.     cudaDeviceSynchronize();
  461.     direct_end = clock();
  463.     //just copy the first result, not all numGpuThreads results
  464.     cudaMemcpy(res, d_res, sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost );
  466.     printf("RESULT: %.16E + %.16E i\n", cuCreal(res[0]), cuCimag(res[0]));
  468.     double timing = (double)((double)(direct_end - direct_start) / CLOCKS_PER_SEC);
  470.     printf("time %.5f s\n", timing);
  472.     cudaDeviceReset();
  474.     return 0;
  475. }

huangapple
  • 本文由 发表于 2023年5月28日 22:03:36
  • 转载请务必保留本文链接：https://go.coder-hub.com/76351879.html
  • cublas
  • cuda
  • determinants
  • eigen
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定