Skip to content

Commit cab0ecd

Browse files
Flying-dragon-boxingFlying-dragon-boxing
authored
Fix CUDA 13+ link error (deepmodeling#7705)
* Update CMake configuration for module_gint * An alternative fix
1 parent c1023ca commit cab0ecd

2 files changed

Lines changed: 227 additions & 326 deletions

File tree

‎source/source_lcao/module_gint/kernel/phi_operator_kernel.cu‎

Lines changed: 9 additions & 321 deletions
Original file line numberDiff line numberDiff line change
@@ -4,234 +4,14 @@
44
#include "source_base/module_device/device.h"
55
#include "source_base/module_device/kernel_compat.h"
66

7-
namespace ModuleGint
8-
{
9-
10-
template<typename Real>
11-
__global__ void set_phi_kernel(
12-
const int nwmax,
13-
const int mgrids_num,
14-
const int nrmax,
15-
const double dr_uniform,
16-
const int* __restrict__ ucell_atom_nwl,
17-
const bool* __restrict__ atom_iw2_new,
18-
const int* __restrict__ atom_iw2_ylm,
19-
const int* __restrict__ atom_nw,
20-
const int* __restrict__ iat2it,
21-
const double* __restrict__ rcut,
22-
const double* __restrict__ psi_u,
23-
const double* __restrict__ dpsi_u,
24-
const double3* __restrict__ mgrids_pos,
25-
const int* __restrict__ atoms_iat,
26-
const double3* __restrict__ atom_rcoords,
27-
const int2* __restrict__ atoms_num_info,
28-
const int* __restrict__ atom_phi_start,
29-
const int* __restrict__ bgrid_phi_len,
30-
Real* __restrict__ phi)
31-
{
32-
const int bgrid_id = blockIdx.y;
33-
const int mgrid_id = blockIdx.x;
34-
const int atoms_num = atoms_num_info[bgrid_id].x;
35-
const int pre_atoms_num = atoms_num_info[bgrid_id].y;
36-
const double3 mgrid_pos = mgrids_pos[mgrid_id];
37-
38-
for (int atom_id = threadIdx.x; atom_id < atoms_num; atom_id += blockDim.x)
39-
{
40-
const int atom_type = iat2it[atoms_iat[atom_id + pre_atoms_num]];
41-
const double3 rcoord = atom_rcoords[atom_id + pre_atoms_num]; // rcoord is the ralative coordinate of an atom and a biggrid
42-
const double3 coord = make_double3(mgrid_pos.x-rcoord.x, // coord is the relative coordinate of an atom and a meshgrid
43-
mgrid_pos.y-rcoord.y,
44-
mgrid_pos.z-rcoord.z);
45-
double dist = norm3d(coord.x, coord.y, coord.z);
46-
if (dist < rcut[atom_type])
47-
{
48-
if (dist < 1.0E-9)
49-
{ dist += 1.0E-9; }
50-
// since nwl is less or equal than 5, the size of ylma is (5+1)^2
51-
double ylma[36];
52-
const int nwl = ucell_atom_nwl[atom_type];
53-
sph_harm(nwl, coord.x/dist, coord.y/dist, coord.z/dist, ylma);
54-
55-
const double pos = dist / dr_uniform;
56-
const int ip = static_cast<int>(pos);
57-
const double dx = pos - ip;
58-
const double dx2 = dx * dx;
59-
const double dx3 = dx2 * dx;
60-
61-
const double c3 = 3.0 * dx2 - 2.0 * dx3;
62-
const double c1 = 1.0 - c3;
63-
const double c2 = (dx - 2.0 * dx2 + dx3) * dr_uniform;
64-
const double c4 = (dx3 - dx2) * dr_uniform;
65-
66-
double psi = 0;
67-
const int it_nw = atom_type * nwmax;
68-
int iw_nr = it_nw * nrmax + ip;
69-
int phi_idx = atom_phi_start[atom_id + pre_atoms_num] +
70-
bgrid_phi_len[bgrid_id] * mgrid_id;
71-
72-
for (int iw = 0; iw < atom_nw[atom_type]; iw++, iw_nr += nrmax)
73-
{
74-
if (atom_iw2_new[it_nw + iw])
75-
{
76-
psi = c1 * psi_u[iw_nr] + c2 * dpsi_u[iw_nr]
77-
+ c3 * psi_u[iw_nr + 1] + c4 * dpsi_u[iw_nr + 1];
78-
}
79-
phi[phi_idx + iw] = static_cast<Real>(psi * ylma[atom_iw2_ylm[it_nw + iw]]);
80-
}
81-
}
82-
else
83-
{
84-
int phi_idx = atom_phi_start[atom_id + pre_atoms_num] +
85-
bgrid_phi_len[bgrid_id] * mgrid_id;
86-
for (int iw = 0; iw < atom_nw[atom_type]; iw++)
87-
{
88-
phi[phi_idx + iw] = Real(0.0);
89-
}
90-
}
91-
}
92-
}
7+
// The template kernels (set_phi_kernel, set_phi_dphi_kernel,
8+
// phi_mul_vldr3_kernel, phi_dot_phi_kernel) are defined in
9+
// phi_operator_kernel.cuh so that their <<<...>>> launches in
10+
// phi_operator_gpu.cu see the definitions in the same translation unit; only
11+
// the non-template kernels live here (see the note in the header).
9312

94-
// Explicit instantiations for set_phi_kernel
95-
template __global__ void set_phi_kernel<double>(
96-
const int, const int, const int, const double,
97-
const int*, const bool*, const int*, const int*, const int*,
98-
const double*, const double*, const double*, const double3*,
99-
const int*, const double3*, const int2*, const int*, const int*,
100-
double*);
101-
template __global__ void set_phi_kernel<float>(
102-
const int, const int, const int, const double,
103-
const int*, const bool*, const int*, const int*, const int*,
104-
const double*, const double*, const double*, const double3*,
105-
const int*, const double3*, const int2*, const int*, const int*,
106-
float*);
107-
108-
template<bool WantPhi>
109-
__global__ void set_phi_dphi_kernel(
110-
const int nwmax,
111-
const int mgrids_num,
112-
const int nrmax,
113-
const double dr_uniform,
114-
const int* __restrict__ ucell_atom_nwl,
115-
const bool* __restrict__ atom_iw2_new,
116-
const int* __restrict__ atom_iw2_ylm,
117-
const int* __restrict__ atom_iw2_l,
118-
const int* __restrict__ atom_nw,
119-
const int* __restrict__ iat2it,
120-
const double* __restrict__ rcut,
121-
const double* __restrict__ psi_u,
122-
const double* __restrict__ dpsi_u,
123-
const double3* __restrict__ mgrids_pos,
124-
const int* __restrict__ atoms_iat,
125-
const double3* __restrict__ atom_rcoords,
126-
const int2* __restrict__ atoms_num_info,
127-
const int* __restrict__ atom_phi_start,
128-
const int* __restrict__ bgrid_phi_len,
129-
double* __restrict__ phi,
130-
double* __restrict__ dphi_x,
131-
double* __restrict__ dphi_y,
132-
double* __restrict__ dphi_z)
13+
namespace ModuleGint
13314
{
134-
const int bgrid_id = blockIdx.y;
135-
const int mgrid_id = blockIdx.x;
136-
const int atoms_num = atoms_num_info[bgrid_id].x;
137-
const int pre_atoms_num = atoms_num_info[bgrid_id].y;
138-
const double3 mgrid_pos = mgrids_pos[mgrid_id];
139-
140-
for (int atom_id = threadIdx.x; atom_id < atoms_num; atom_id += blockDim.x)
141-
{
142-
const int atom_type = iat2it[atoms_iat[atom_id + pre_atoms_num]];
143-
const double3 rcoord = atom_rcoords[atom_id + pre_atoms_num];
144-
const double3 coord = make_double3(mgrid_pos.x-rcoord.x,
145-
mgrid_pos.y-rcoord.y,
146-
mgrid_pos.z-rcoord.z);
147-
double dist = norm3d(coord.x, coord.y, coord.z);
148-
if (dist < rcut[atom_type])
149-
{
150-
if (dist < 1.0E-9)
151-
{ dist += 1.0E-9; }
152-
// since nwl is less or equal than 5, the size of rly is (5+1)^2
153-
// size of grly = 36 * 3
154-
double rly[36];
155-
double grly[36 * 3];
156-
const int nwl = ucell_atom_nwl[atom_type];
157-
grad_rl_sph_harm(nwl, coord.x, coord.y, coord.z, rly, grly);
158-
159-
// interpolation
160-
const double inv_dist = 1.0 / dist; // hoisted: re-used by every iw below
161-
const double pos = dist / dr_uniform;
162-
const int ip = static_cast<int>(pos);
163-
const double x0 = pos - ip;
164-
const double x1 = 1.0 - x0;
165-
const double x2 = 2.0 - x0;
166-
const double x3 = 3.0 - x0;
167-
const double x12 = x1 * x2 / 6;
168-
const double x03 = x0 * x3 / 2;
169-
double tmp = 0;
170-
double dtmp = 0;
171-
const int it_nw = atom_type * nwmax;
172-
int iw_nr = it_nw * nrmax + ip;
173-
int phi_idx = atom_phi_start[atom_id + pre_atoms_num] +
174-
bgrid_phi_len[bgrid_id] * mgrid_id;
175-
for (int iw = 0; iw < atom_nw[atom_type]; iw++, iw_nr += nrmax)
176-
{
177-
if (atom_iw2_new[it_nw + iw])
178-
{
179-
tmp = x12 * (psi_u[iw_nr] * x3 + psi_u[iw_nr + 3] * x0)
180-
+ x03 * (psi_u[iw_nr + 1] * x2 - psi_u[iw_nr + 2] * x1);
181-
dtmp = x12 * (dpsi_u[iw_nr] * x3 + dpsi_u[iw_nr + 3] * x0)
182-
+ x03 * (dpsi_u[iw_nr + 1] * x2 - dpsi_u[iw_nr + 2] * x1);
183-
}
184-
const int iw_l = atom_iw2_l[it_nw + iw];
185-
const int idx_ylm = atom_iw2_ylm [it_nw + iw];
186-
const double rl = pow_int(dist, iw_l);
187-
const double inv_rl = 1.0 / rl;
188-
const double tmprl = tmp * inv_rl;
189-
190-
if (WantPhi)
191-
{
192-
phi[phi_idx + iw] = tmprl * rly[idx_ylm];
193-
}
194-
// derivative of wave functions with respect to atom positions.
195-
// (dtmp - tmp*iw_l/dist) / rl * rly / dist == (dtmp*inv_dist - tmp*iw_l*inv_dist^2) * inv_rl * rly
196-
const double tmpdphi_rly = (dtmp * inv_dist - tmp * iw_l * inv_dist * inv_dist)
197-
* inv_rl * rly[idx_ylm];
198-
199-
dphi_x[phi_idx + iw] = tmpdphi_rly * coord.x + tmprl * grly[idx_ylm * 3 + 0];
200-
dphi_y[phi_idx + iw] = tmpdphi_rly * coord.y + tmprl * grly[idx_ylm * 3 + 1];
201-
dphi_z[phi_idx + iw] = tmpdphi_rly * coord.z + tmprl * grly[idx_ylm * 3 + 2];
202-
}
203-
}
204-
else
205-
{
206-
int phi_idx = atom_phi_start[atom_id + pre_atoms_num] +
207-
bgrid_phi_len[bgrid_id] * mgrid_id;
208-
for (int iw = 0; iw < atom_nw[atom_type]; iw++)
209-
{
210-
if (WantPhi)
211-
{
212-
phi[phi_idx + iw] = 0.0;
213-
}
214-
dphi_x[phi_idx + iw] = 0.0;
215-
dphi_y[phi_idx + iw] = 0.0;
216-
dphi_z[phi_idx + iw] = 0.0;
217-
}
218-
}
219-
}
220-
}
221-
222-
// Explicit instantiations for set_phi_dphi_kernel
223-
template __global__ void set_phi_dphi_kernel<true>(
224-
const int, const int, const int, const double,
225-
const int*, const bool*, const int*, const int*, const int*, const int*,
226-
const double*, const double*, const double*, const double3*,
227-
const int*, const double3*, const int2*, const int*, const int*,
228-
double*, double*, double*, double*);
229-
template __global__ void set_phi_dphi_kernel<false>(
230-
const int, const int, const int, const double,
231-
const int*, const bool*, const int*, const int*, const int*, const int*,
232-
const double*, const double*, const double*, const double3*,
233-
const int*, const double3*, const int2*, const int*, const int*,
234-
double*, double*, double*, double*);
23515

23616
// The code for `set_ddphi_kernel` is quite difficult to understand.
23717
// To grasp it, you better refer to the CPU function `set_ddphi`
@@ -267,7 +47,7 @@ __global__ void set_ddphi_kernel(
26747
const int atoms_num = atoms_num_info[bgrid_id].x;
26848
const int pre_atoms_num = atoms_num_info[bgrid_id].y;
26949
const double3 mgrid_pos = mgrids_pos[mgrid_id];
270-
50+
27151
for (int atom_id = threadIdx.x; atom_id < atoms_num; atom_id += blockDim.x)
27252
{
27353
const int atom_type = iat2it[atoms_iat[atom_id + pre_atoms_num]];
@@ -324,7 +104,7 @@ __global__ void set_ddphi_kernel(
324104
const double tmprl = tmp * inv_rl;
325105
const double tmpdphi_rly = (dtmp * inv_dist - tmp * iw_l * inv_dist * inv_dist)
326106
* inv_rl * rly[idx_ylm];
327-
107+
328108
double dphi[3];
329109
dphi[0] = tmpdphi_rly * coord[0] + tmprl * grly[idx_ylm * 3 + 0];
330110
dphi[1] = tmpdphi_rly * coord[1] + tmprl * grly[idx_ylm * 3 + 1];
@@ -378,98 +158,6 @@ __global__ void set_ddphi_kernel(
378158
}
379159
}
380160

381-
template<typename Real>
382-
__global__ void phi_mul_vldr3_kernel(
383-
const Real* __restrict__ vl,
384-
const Real dr3,
385-
const Real* __restrict__ phi,
386-
const int mgrids_per_bgrid,
387-
const int* __restrict__ mgrid_lidx,
388-
const int* __restrict__ bgrid_phi_len,
389-
const int* __restrict__ bgrid_phi_start,
390-
Real* __restrict__ result)
391-
{
392-
const int bgrid_id = blockIdx.y;
393-
const int mgrid_id = blockIdx.x;
394-
const int phi_len = bgrid_phi_len[bgrid_id];
395-
const int phi_start = bgrid_phi_start[bgrid_id] + mgrid_id * phi_len;
396-
const int batch_mgrid_id = bgrid_id * mgrids_per_bgrid + mgrid_id;
397-
const Real vldr3 = vl[mgrid_lidx[batch_mgrid_id]] * dr3;
398-
for(int i = threadIdx.x; i < phi_len; i += blockDim.x)
399-
{
400-
result[phi_start + i] = phi[phi_start + i] * vldr3;
401-
}
402-
}
403-
404-
// Explicit instantiations for phi_mul_vldr3_kernel
405-
template __global__ void phi_mul_vldr3_kernel<double>(
406-
const double*, const double, const double*, const int,
407-
const int*, const int*, const int*, double*);
408-
template __global__ void phi_mul_vldr3_kernel<float>(
409-
const float*, const float, const float*, const int,
410-
const int*, const int*, const int*, float*);
411-
412-
// rho(ir) = \sum_{iwt} \phi_i(ir,iwt) * \phi_j^*(ir,iwt)
413-
// each block calculate the dot product of phi_i and phi_j of a meshgrid.
414-
// The per-thread/warp/block reduction is in double regardless of input types
415-
// so that fp32 inputs are summed without catastrophic precision loss.
416-
template<typename Tin_a, typename Tin_b>
417-
__global__ void phi_dot_phi_kernel(
418-
const Tin_a* __restrict__ phi_i,
419-
const Tin_b* __restrict__ phi_j,
420-
const int mgrids_per_bgrid,
421-
const int* __restrict__ mgrid_lidx,
422-
const int* __restrict__ bgrid_phi_len,
423-
const int* __restrict__ bgrid_phi_start,
424-
double* __restrict__ rho)
425-
{
426-
__shared__ double s_data[32]; // the length of s_data equals the max warp num of a block
427-
const int bgrid_id = blockIdx.y;
428-
const int mgrid_id = blockIdx.x;
429-
const int phi_len = bgrid_phi_len[bgrid_id];
430-
const int phi_start = bgrid_phi_start[bgrid_id] + mgrid_id * phi_len;
431-
const Tin_a* phi_i_mgrid = phi_i + phi_start;
432-
const Tin_b* phi_j_mgrid = phi_j + phi_start;
433-
const int batch_mgrid_id = bgrid_id * mgrids_per_bgrid + mgrid_id;
434-
const int mgrid_local_idx = mgrid_lidx[batch_mgrid_id];
435-
const int tid = threadIdx.x;
436-
const int warp_id = tid / 32;
437-
const int lane_id = tid % 32;
438-
double tmp_sum = 0.0;
439-
440-
for (int i = tid; i < phi_len; i += blockDim.x)
441-
{
442-
tmp_sum += phi_i_mgrid[i] * phi_j_mgrid[i];
443-
}
444-
445-
tmp_sum = warpReduceSum(tmp_sum);
446-
447-
if (lane_id == 0)
448-
{
449-
s_data[warp_id] = tmp_sum;
450-
}
451-
__syncthreads();
452-
453-
tmp_sum = (tid < blockDim.x / 32) ? s_data[tid] : 0.0;
454-
if(warp_id == 0)
455-
{
456-
tmp_sum = warpReduceSum(tmp_sum);
457-
}
458-
459-
if(tid == 0)
460-
{
461-
atomicAdd(&rho[mgrid_local_idx], tmp_sum);
462-
}
463-
}
464-
465-
// Explicit instantiations for phi_dot_phi_kernel
466-
template __global__ void phi_dot_phi_kernel<double, double>(
467-
const double*, const double*, const int,
468-
const int*, const int*, const int*, double*);
469-
template __global__ void phi_dot_phi_kernel<float, double>(
470-
const float*, const double*, const int,
471-
const int*, const int*, const int*, double*);
472-
473161
__global__ void phi_dot_dphi_kernel(
474162
const double* __restrict__ phi,
475163
const double* __restrict__ dphi_x,
@@ -580,7 +268,7 @@ __global__ void phi_dot_dphi_r_kernel(
580268
}
581269
}
582270
}
583-
271+
584272
// single-warp reduce
585273
#pragma unroll
586274
for (int i = 0; i < 6; i++)

0 commit comments

Comments
 (0)