| 92 | |
| 93 | |
| 94 | __global__ void DeviceSparseCholeskySolveLKernel ( |
| 95 | FlatTable<int> dependency, |
| 96 | FlatArray<Dev<int>> incomingdep, |
| 97 | FlatVector<Dev<double>> hy, |
| 98 | int & cnt, |
| 99 | FlatArray<Dev<MicroTask>> microtasks, |
| 100 | FlatArray<Dev<int>> blocks, |
| 101 | FlatArray<Dev<int>> rowindex2, |
| 102 | FlatArray<Dev<size_t>> firstinrow_ri, |
| 103 | FlatArray<Dev<size_t>> firstinrow, |
| 104 | FlatArray<Dev<double>> lfact |
| 105 | ) |
| 106 | { |
| 107 | __shared__ int myjobs[16]; // max blockDim.y; |
| 108 | DeviceBlockRegionTracer brt(gridDim.x*blockDim.y, gridDim.x*threadIdx.y + blockIdx.x, threadIdx.x); |
| 109 | |
| 110 | while (true) |
| 111 | { |
| 112 | if (threadIdx.x == 0) |
| 113 | myjobs[threadIdx.y] = atomicAdd(&cnt, 1); |
| 114 | __syncwarp(); |
| 115 | |
| 116 | int myjob = myjobs[threadIdx.y]; |
| 117 | if (myjob >= dependency.Size()) |
| 118 | break; |
| 119 | |
| 120 | volatile int * n_deps = (int*)&incomingdep[myjob]; |
| 121 | if (threadIdx.x == 0) |
| 122 | while(*n_deps); |
| 123 | __syncwarp(); |
| 124 | |
| 125 | // do the work for myjob.... |
| 126 | |
| 127 | MicroTask task = microtasks[myjob]; |
| 128 | size_t blocknr = task.blocknr; |
| 129 | auto range = T_Range<int>(blocks[blocknr], blocks[blocknr+1]); |
| 130 | auto base = firstinrow_ri[range.First()] + range.Size()-1; |
| 131 | auto ext_size = firstinrow[range.First()+1]-firstinrow[range.First()] - range.Size()+1; |
| 132 | auto extdofs = rowindex2.Range(base, base+ext_size); |
| 133 | |
| 134 | __threadfence(); |
| 135 | |
| 136 | if ((task.type == MicroTask::L_BLOCK) || (task.type == MicroTask::LB_BLOCK)) |
| 137 | { |
| 138 | DeviceRegionTracer rt(brt, 0, task.blocknr); |
| 139 | for (int i = blocks[blocknr]; i < blocks[blocknr+1]-1; i++) |
| 140 | { |
| 141 | size_t size = range.end()-i-1; |
| 142 | |
| 143 | auto vlfact = lfact.Range(firstinrow[i]-i-1, firstinrow[i]+size); |
| 144 | |
| 145 | __threadfence_block(); |
| 146 | |
| 147 | double hyi = hy(i); |
| 148 | for (int j = threadIdx.x+blocks[blocknr]; j < range.end(); j += blockDim.x) |
| 149 | if (j > i) |
| 150 | hy(j) -= vlfact[j] * hyi; |
| 151 | __threadfence_block(); |