| 50 | |
| 51 | template<typename T> |
| 52 | __device__ void JacobiSVD(int m, int n) { |
| 53 | const int iterations = 30; |
| 54 | |
| 55 | int tid_x = threadIdx.x; |
| 56 | int bsz_x = blockDim.x; |
| 57 | int tid_y = threadIdx.y; |
| 58 | // int gid_y = blockIdx.y * blockDim.y + tid_y; |
| 59 | |
| 60 | __shared__ T s_acc1[256]; |
| 61 | __shared__ T s_acc2[256]; |
| 62 | |
| 63 | __shared__ T s_d[16 * 9]; |
| 64 | |
| 65 | T* s_V = (T*)sh; |
| 66 | T* s_S = (T*)sh + 16 * 81; |
| 67 | |
| 68 | int doff = tid_y * n; |
| 69 | int soff = tid_y * 81; |
| 70 | |
| 71 | if (tid_x < n) { |
| 72 | T acc1 = 0; |
| 73 | for (int i = 0; i < m; i++) { |
| 74 | int stid = soff + tid_x * m + i; |
| 75 | T t = s_S[stid]; |
| 76 | acc1 += t * t; |
| 77 | s_V[stid] = (tid_x == i) ? 1 : 0; |
| 78 | } |
| 79 | s_d[doff + tid_x] = acc1; |
| 80 | } |
| 81 | __syncthreads(); |
| 82 | |
| 83 | for (int it = 0; it < iterations; it++) { |
| 84 | for (int i = 0; i < n - 1; i++) { |
| 85 | for (int j = i + 1; j < n; j++) { |
| 86 | T* Si = s_S + soff + i * m; |
| 87 | T* Sj = s_S + soff + j * m; |
| 88 | |
| 89 | T* Vi = s_V + soff + i * n; |
| 90 | T* Vj = s_V + soff + j * n; |
| 91 | |
| 92 | T p = (T)0; |
| 93 | for (int k = 0; k < m; k++) p += Si[k] * Sj[k]; |
| 94 | |
| 95 | T di = s_d[doff + i]; |
| 96 | T dj = s_d[doff + j]; |
| 97 | __syncthreads(); |
| 98 | |
| 99 | T c = 0, s = 0; |
| 100 | T t0 = 0, t1 = 0; |
| 101 | int cond = (fabs(p) > m * EPS<T>::eps() * sqrt(di * dj)); |
| 102 | T a = 0, b = 0; |
| 103 | |
| 104 | if (cond) { |
| 105 | T y = di - dj; |
| 106 | T r = hypot(p * 2, y); |
| 107 | T r2 = r * 2; |
| 108 | if (y >= 0) { |
| 109 | c = sqrt((r + y) / r2); |