Skip to content

Commit b1c36c1

Browse files
authored
Batch DF derivative work and reduce gradient workspace (#967)
* df/grad/td: stabilize metric application with factor products Improve numerical stability of the density-pair gradient against changes in memory-dependent contraction sizes. For density pairs from test_jk_energy_per_atom_dm_pairs in gpu4pyscf/df/tests/test_df_tdrhf_grad.py, substituting 16.2 MB and 2.49 MB as get_avail_mem returns produces a maximum K gradient component discrepancy of 4e-12 Eh/bohr with the expanded metric. _jk_energies_by_dm_factors forms G = C M^-1 C.T and applies G X to assembled J vectors and K tensors. M is the contracted auxiliary Coulomb matrix; C expands its basis into the integral representation. Small eigenvalues of M produce large entries in G, so its output dot products can subtract large terms to obtain small coefficients. Compute t = A.T X, then A t, for each assembled array X, sharing A across J and K. Following gpu4pyscf.df.df._decompose_j2c, use A = C L^-T for Cholesky M = L L.T, or A = C U lambda^-1/2 with retained eigenvectors U and eigenvalues lambda. Preserve decomposition selection, fallback, and the eigenvalue cutoff. In this fixture, G entries reach about 20000, while A entries reach about 70. One coefficient near -0.0065 sums positive and negative contributions of about 50000 in G X. Along the factorized path, a contributing t = A.T X entry combines signed sums of about 200 to give an intermediate near -0.3; the final A t combines signed sums of about 20 to give the coefficient near -0.0065. The first product thus performs part of the cancellation before forming the output coefficients. The final product reuses these smaller t values across outputs. An error in one t value enters each output with the corresponding weight from A; when those contributions cancel in the gradient, their shared error cancels too. Dense G X computes a separate large-term sum for each output, with its own rounding residual. Both factorized products still round; the gradient comparison below measures their combined effect. For the same memory probes, the factorized application reduces the maximum K gradient component discrepancy to 2e-14 Eh/bohr. Multiplying the same stored G in extended precision also reduces the discrepancy to 1e-14, isolating rounding in the dense product. The improved consistency trades one dense metric multiplication for two factor multiplications. With n expanded and m contracted auxiliary functions, Cholesky application uses 2*m/n times the leading arithmetic of G X; eigenvalue truncation replaces m with the retained rank. * df/grad: size scratch for auxiliary tiles and density factors Memory-derived auxiliary sizes can exceed the evaluator's largest tile. Cap integral and contraction scratch at that tile size in RHF, UHF, and the TD summed-density and density-pair paths. The same cap also bounds derivative batching. Hermitian exchange scratch also reserves excess AO-sized capacity. Size it by the density-factor width and AO count in RHF, UHF, and summed-density TD, retaining both spin channels in UHF. Preserve the transpose workspace required by non-Hermitian inputs. * df/grad,hessian: pack auxiliary groups in derivative launches RHF, UHF and TD-RHF derivative drivers filled one auxiliary group per launch despite buffers sized for several groups, limiting grid parallelism. Pack consecutive groups within the existing auxiliary-function capacity and pass the group count through the native driver's existing second grid axis. This exposes more independent blocks per packed launch while preserving group boundaries and function ordering. To keep RHF Hessian packed work distributable across devices, target the auxiliary-function count per device when limiting tile size, with whole groups setting the minimum granularity. Assign packed tiles through the existing shared iterator. * df/int3c2e_bdiv: pack auxiliary batches Int3c2eOpt.int3c2e_evaluator evaluates three-center Coulomb integrals. Its aux_batch_size argument specifies a requested auxiliary-function capacity, but setting it produced one shell group per call even when several groups fit. This left buffer capacity and GPU parallelism unused. Use splits_by_blocksize to pack consecutive groups using that capacity. Apply the same rule to _int3c2e_ip1_evaluator, which evaluates first derivatives of the three-center integrals, and pass the packed group count to its native driver. Preserve auxiliary-function ordering; both evaluators return offsets describing their packed batches. Larger tiles expose an existing buffer-sizing bug in tdscf.ris.get_Tpq. It allocated the integral-output buffer using the smaller inner orbital-contraction block size, causing buffer assertions under limited memory. Fix this by allocating for the largest returned integral tile, then contracting it in smaller orbital blocks. Add test_get_Tpq_limited_memory in df/tests/test_df_tddft_ris.py. It compares J/K orbital transforms under normal and constrained memory, covering both host and GPU result storage, to guard the full-tile buffer allocation.
1 parent 9af748a commit b1c36c1

9 files changed

Lines changed: 118 additions & 47 deletions

File tree

‎gpu4pyscf/df/grad/rhf.py‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
from cupyx.scipy.linalg import solve_triangular
1919
from pyscf import lib
2020
from gpu4pyscf.lib import logger
21+
from gpu4pyscf.lib.utils import splits_by_blocksize
2122
from gpu4pyscf.lib.cupy_helper import (
2223
contract, asarray, ndarray, cholesky, eigh, transpose_sum, get_avail_mem)
2324
from gpu4pyscf.grad import rhf as rhf_grad
@@ -98,9 +99,10 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
9899
eval_j3c, aux_sorting, _, aux_offsets = int3c2e_opt.int3c2e_evaluator(
99100
aux_batch_size=batch_size, reorder_aux=True, cart=True,
100101
omega=omega, lr_factor=lr_factor, sr_factor=sr_factor)
102+
batch_size = min(batch_size, int((aux_offsets[1:]-aux_offsets[:-1]).max()))
101103
aux_batches = len(aux_offsets) - 1
102104

103-
blksize = max(1, min(naux, int(mem_avail*.4/(nao*nao*2*8))//8*8))
105+
blksize = max(1, min(batch_size, int(mem_avail*.4/(nao*nao*2*8))//8*8))
104106
log.debug1('%.3f GB free memory. nao_pair=%d naux=%d batch_size=%d blksize=%d',
105107
mem_free*1e-9, nao_pair, naux, batch_size, blksize)
106108

@@ -169,11 +171,13 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
169171
aux0 = aux1 = 0
170172
buf = cp.empty((nao_pair*batch_size))
171173
buf1 = cp.empty((blksize, nao, nao))
172-
buf2 = cp.empty((blksize, nao, nao))
174+
# Only non-Hermitian densities reuse this buffer for an AO-sized transpose.
175+
buf2 = cp.empty((blksize, nocc if hermi == 1 else nao, nao))
173176
ejk = cp.zeros((mol.natm, 3))
174177
ejk_aux = cp.zeros((mol.natm, 3))
175-
for kbatch in range(aux_batches):
176-
naux_in_batch = aux_offsets[kbatch+1] - aux_offsets[kbatch]
178+
aux_splits = splits_by_blocksize(aux_offsets, batch_size)
179+
for kbatch, kend in zip(aux_splits[:-1], aux_splits[1:]):
180+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
177181
aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
178182
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf)
179183
for k0, k1 in lib.prange(0, naux_in_batch, blksize):
@@ -208,7 +212,7 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
208212
ctypes.byref(int3c2e_envs),
209213
ctypes.c_int(shm_size_max),
210214
ctypes.c_int(len(shl_pair_offsets) - 1),
211-
ctypes.c_int(1),
215+
ctypes.c_int(kend-kbatch),
212216
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
213217
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
214218
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),

‎gpu4pyscf/df/grad/tdrhf.py‎

Lines changed: 43 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,12 @@
1616
import ctypes
1717
import numpy as np
1818
import cupy as cp
19+
from cupyx.scipy.linalg import solve_triangular
1920
from pyscf import lib
2021
from gpu4pyscf.lib import logger
22+
from gpu4pyscf.lib.utils import splits_by_blocksize
2123
from gpu4pyscf.lib.cupy_helper import (
22-
contract, asarray, ndarray, transpose_sum, get_avail_mem)
24+
contract, asarray, ndarray, cholesky, eigh, transpose_sum, get_avail_mem)
2325
from gpu4pyscf.df.grad.rhf import (
2426
_split_l_ctr_pattern, get_ao_pair_loc, libvhf_rys, Int3c2eOpt, int2c2e,
2527
int3c2e_scheme, _gen_metric_solver, _factorize_dm)
@@ -75,9 +77,10 @@ def _jk_energy_per_atom(int3c2e_opt, dms, j_factor=None, k_factor=None, hermi=0,
7577
eval_j3c, aux_sorting, _, aux_offsets = int3c2e_opt.int3c2e_evaluator(
7678
aux_batch_size=batch_size, reorder_aux=True, cart=True,
7779
omega=omega, lr_factor=lr_factor, sr_factor=sr_factor)
80+
batch_size = min(batch_size, int((aux_offsets[1:]-aux_offsets[:-1]).max()))
7881
aux_batches = len(aux_offsets) - 1
7982

80-
blksize = max(1, min(naux, int(mem_avail*.4/(nao*nao*2*8))//8*8))
83+
blksize = max(1, min(batch_size, int(mem_avail*.4/(nao*nao*2*8))//8*8))
8184
log.debug('%.3f GB free memory. nao_pair=%d naux=%d batch_size=%d blksize=%d',
8285
mem_free*1e-9, nao_pair, naux, batch_size, blksize)
8386

@@ -142,23 +145,23 @@ def _jk_energy_per_atom(int3c2e_opt, dms, j_factor=None, k_factor=None, hermi=0,
142145
l_ctr_aux_offsets, auxmol.uniq_l_ctr, batch_size)
143146
ksh_offsets_cpu = l_ctr_aux_offsets
144147
ksh_offsets_gpu = cp.asarray(ksh_offsets_cpu+mol.nbas, dtype=np.int32)
145-
l_ctr_aux_counts = l_ctr_aux_offsets[1:] - l_ctr_aux_offsets[:-1]
148+
aux_offsets = aux_loc[ksh_offsets_cpu]
146149

147150
if j_factor is not None:
148151
dms = contract('npi,nqi->npq', dm_factor_l, dm_factor_r)
149152

150153
int3c2e_envs = int3c2e_opt.int3c2e_envs
151154
kern = libvhf_rys.sum_ejk_int3c2e_ip1
152-
l = np.arange(laux+1)
153-
nf = (l + 1) * (l + 2) // 2
154155
aux0 = aux1 = 0
155156
buf = cp.empty((nao_pair*batch_size))
156157
buf1 = cp.empty((blksize, nao, nao))
157-
buf2 = cp.empty((blksize, nao, nao))
158+
# Non-Hermitian densities also use this buffer for an AO-sized transpose.
159+
buf2 = cp.empty((blksize, nocc if hermi == 1 else nao, nao))
158160
ejk = cp.zeros((mol.natm, 3))
159161
ejk_aux = cp.zeros((mol.natm, 3))
160-
for kbatch, lk, in enumerate(uniq_l_ctr_aux[:,0]):
161-
naux_in_batch = nf[lk] * l_ctr_aux_counts[kbatch]
162+
aux_splits = splits_by_blocksize(aux_offsets, batch_size)
163+
for kbatch, kend in zip(aux_splits[:-1], aux_splits[1:]):
164+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
162165
aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
163166
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf)
164167
for k0, k1 in lib.prange(0, naux_in_batch, blksize):
@@ -194,7 +197,7 @@ def _jk_energy_per_atom(int3c2e_opt, dms, j_factor=None, k_factor=None, hermi=0,
194197
ctypes.byref(int3c2e_envs),
195198
ctypes.c_int(shm_size_max),
196199
ctypes.c_int(len(shl_pair_offsets) - 1),
197-
ctypes.c_int(1),
200+
ctypes.c_int(kend-kbatch),
198201
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
199202
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
200203
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),
@@ -441,9 +444,10 @@ def _jk_energies_by_dm_factors(int3c2e_opt, dm_factors, j_factor, k_factor,
441444
eval_j3c, aux_sorting, _, aux_offsets = int3c2e_opt.int3c2e_evaluator(
442445
aux_batch_size=batch_size, reorder_aux=True, cart=True,
443446
omega=omega, lr_factor=lr_factor, sr_factor=sr_factor)
447+
batch_size = min(batch_size, int((aux_offsets[1:]-aux_offsets[:-1]).max()))
444448
aux_batches = len(aux_offsets) - 1
445449

446-
blksize = max(1, min(naux, int(mem_avail*.45/(nao*nao*2*8))//8*8))
450+
blksize = max(1, min(batch_size, int(mem_avail*.45/(nao*nao*2*8))//8*8))
447451
log.debug('%.3f GB free memory. nao_pair=%d naux=%d batch_size=%d blksize=%d',
448452
mem_free*1e-9, nao_pair, naux, batch_size, blksize)
449453

@@ -481,23 +485,39 @@ def _jk_energies_by_dm_factors(int3c2e_opt, dm_factors, j_factor, k_factor,
481485
tmp = None
482486

483487
j2c = int2c2e(auxmol, omega=omega, lr_factor=lr_factor, sr_factor=sr_factor)
488+
# Factor C M^-1 C.T as A A.T, following df.df._decompose_j2c,
489+
# while retaining the derivative builder's explicit CD/ED selection.
484490
if mol.omega <= 0 and not auxmol.mol.cart:
485-
metric = aux_coeff.dot(_gen_metric_solver(j2c, 'CD')(aux_coeff.T))
491+
decompose_j2c = 'CD'
492+
try:
493+
cd_low = cholesky(j2c)
494+
aux_coeff = solve_triangular(
495+
cd_low, aux_coeff.T, lower=True, overwrite_b=True).T
496+
cd_low = None
497+
except RuntimeError:
498+
decompose_j2c = 'ED'
486499
else:
487-
metric = aux_coeff.dot(_gen_metric_solver(j2c, 'ED')(aux_coeff.T))
488-
j2c = aux_coeff = None
500+
decompose_j2c = 'ED'
501+
if decompose_j2c == 'ED':
502+
w, v = eigh(j2c)
503+
idx = cp.where(w > df.LINEAR_DEP_THR)[0]
504+
aux_coeff = aux_coeff.dot(v[:,idx] / cp.sqrt(w[idx]))
505+
w = v = idx = None
489506
for i in range(n_dm):
490-
j3c_o2o1[i] = contract('uv,vij->uij', metric, j3c_o2o1[i])
491-
j3c_o1o2[i] = contract('uv,vij->uij', metric, j3c_o1o2[i])
507+
tmp = contract('up,uij->pij', aux_coeff, j3c_o2o1[i])
508+
j3c_o2o1[i] = contract('up,pij->uij', aux_coeff, tmp)
509+
tmp = contract('up,uij->pij', aux_coeff, j3c_o1o2[i])
510+
j3c_o1o2[i] = contract('up,pij->uij', aux_coeff, tmp)
511+
j2c = tmp = None
492512
cp.get_default_memory_pool().free_all_blocks()
493513

494514
if j_factor is not None:
495515
j_factor = cp.asarray(j_factor)
496-
auxvec1 = cp.einsum('uv,nv->nu', metric, auxvec1)
497-
auxvec2 = cp.einsum('uv,nv->nu', metric, auxvec2)
516+
auxvec1 = aux_coeff.dot(aux_coeff.T.dot(auxvec1.T)).T
517+
auxvec2 = aux_coeff.dot(aux_coeff.T.dot(auxvec2.T)).T
498518
auxvec1_jfac = j_factor[:,None] * auxvec1
499519
auxvec2_jfac = j_factor[:,None] * auxvec2
500-
metric = None
520+
aux_coeff = None
501521

502522
# (d/dX P|Q) contributions
503523
if sum_results:
@@ -552,11 +572,9 @@ def _jk_energies_by_dm_factors(int3c2e_opt, dm_factors, j_factor, k_factor,
552572
l_ctr_aux_offsets, auxmol.uniq_l_ctr, batch_size)
553573
ksh_offsets_cpu = l_ctr_aux_offsets
554574
ksh_offsets_gpu = cp.asarray(ksh_offsets_cpu+mol.nbas, dtype=np.int32)
555-
l_ctr_aux_counts = l_ctr_aux_offsets[1:] - l_ctr_aux_offsets[:-1]
575+
aux_offsets = aux_loc[ksh_offsets_cpu]
556576

557577
int3c2e_envs = int3c2e_opt.int3c2e_envs
558-
l = np.arange(laux+1)
559-
nf = (l + 1) * (l + 2) // 2
560578
aux0 = aux1 = 0
561579
buf1 = cp.empty((blksize, nao, nao))
562580
buf2 = cp.empty((blksize, nao, nao))
@@ -569,8 +587,9 @@ def _jk_energies_by_dm_factors(int3c2e_opt, dm_factors, j_factor, k_factor,
569587
ejk = cp.zeros((n_dm, mol.natm, 3))
570588
buf = cp.empty((n_dm*nao_pair*batch_size))
571589

572-
for kbatch, lk, in enumerate(uniq_l_ctr_aux[:,0]):
573-
naux_in_batch = nf[lk] * l_ctr_aux_counts[kbatch]
590+
aux_splits = splits_by_blocksize(aux_offsets, batch_size)
591+
for kbatch, kend in zip(aux_splits[:-1], aux_splits[1:]):
592+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
574593
aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
575594
if sum_results:
576595
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf3)
@@ -616,7 +635,7 @@ def _jk_energies_by_dm_factors(int3c2e_opt, dm_factors, j_factor, k_factor,
616635
ctypes.byref(int3c2e_envs),
617636
ctypes.c_int(shm_size_max),
618637
ctypes.c_int(len(shl_pair_offsets) - 1),
619-
ctypes.c_int(1),
638+
ctypes.c_int(kend-kbatch),
620639
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
621640
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
622641
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),

‎gpu4pyscf/df/grad/uhf.py‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
import cupy as cp
1818
from pyscf import lib
1919
from gpu4pyscf.lib import logger
20+
from gpu4pyscf.lib.utils import splits_by_blocksize
2021
from gpu4pyscf.lib.cupy_helper import contract, asarray, ndarray, get_avail_mem
2122
from gpu4pyscf.grad import uhf as uhf_grad
2223
from gpu4pyscf.df.grad.rhf import (
@@ -73,9 +74,10 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
7374
eval_j3c, aux_sorting, _, aux_offsets = int3c2e_opt.int3c2e_evaluator(
7475
aux_batch_size=batch_size, reorder_aux=True, cart=True,
7576
omega=omega, lr_factor=lr_factor, sr_factor=sr_factor)
77+
batch_size = min(batch_size, int((aux_offsets[1:]-aux_offsets[:-1]).max()))
7678
aux_batches = len(aux_offsets) - 1
7779

78-
blksize = max(1, min(naux, int(mem_avail*.4/(nao*nao*2*8))//8*8))
80+
blksize = max(1, min(batch_size, int(mem_avail*.4/(nao*nao*2*8))//8*8))
7981
log.debug1('%.3f GB free memory. nao_pair=%d naux=%d batch_size=%d blksize=%d',
8082
mem_free*1e-9, nao_pair, naux, batch_size, blksize)
8183

@@ -145,11 +147,13 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
145147
aux0 = aux1 = 0
146148
buf = cp.empty((nao_pair*batch_size))
147149
buf1 = cp.empty((blksize, nao, nao))
148-
buf2 = cp.empty(blksize*nao*max(2*nocc, nao))
150+
# Non-Hermitian densities also use this buffer for an AO-sized transpose.
151+
buf2 = cp.empty(blksize*nao*(2*nocc if hermi == 1 else max(2*nocc, nao)))
149152
ejk = cp.zeros((mol.natm, 3))
150153
ejk_aux = cp.zeros((mol.natm, 3))
151-
for kbatch in range(aux_batches):
152-
naux_in_batch = aux_offsets[kbatch+1] - aux_offsets[kbatch]
154+
aux_splits = splits_by_blocksize(aux_offsets, batch_size)
155+
for kbatch, kend in zip(aux_splits[:-1], aux_splits[1:]):
156+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
153157
aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
154158
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf)
155159
for k0, k1 in lib.prange(0, naux_in_batch, blksize):
@@ -184,7 +188,7 @@ def _jk_energy_per_atom(int3c2e_opt, dm, j_factor=1, k_factor=1, hermi=0,
184188
ctypes.byref(int3c2e_envs),
185189
ctypes.c_int(shm_size_max),
186190
ctypes.c_int(len(shl_pair_offsets) - 1),
187-
ctypes.c_int(1),
191+
ctypes.c_int(kend-kbatch),
188192
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
189193
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
190194
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),

‎gpu4pyscf/df/hessian/rhf.py‎

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
from pyscf import lib
2828
from pyscf.gto import ATOM_OF
2929
from gpu4pyscf.lib import logger
30+
from gpu4pyscf.lib.utils import splits_by_blocksize
3031
from gpu4pyscf.lib.cupy_helper import (
3132
contract, ndarray, transpose_sum, get_avail_mem, empty_mapped, fill_symmetric)
3233
from gpu4pyscf.df.int3c2e_bdiv import (
@@ -199,7 +200,12 @@ def proc(aux_batch_iter):
199200
ksh_offsets_cpu = l_ctr_aux_offsets
200201
ksh_offsets_gpu = cp.asarray(ksh_offsets_cpu+mol.nbas, dtype=np.int32)
201202
aux_offsets = aux_loc[ksh_offsets_cpu]
202-
num_aux_batches = len(aux_offsets) - 1
203+
# Keep device tasks small enough to distribute, retaining the existing
204+
# group boundaries that define the auxiliary function ordering.
205+
aux_batch_size = min(batch_size, max(np.diff(aux_offsets).max(),
206+
(naux+num_devices-1)//num_devices))
207+
aux_splits = splits_by_blocksize(aux_offsets, aux_batch_size)
208+
num_aux_batches = len(aux_splits) - 1
203209

204210
# (20|0)(0|0)(0|00) + (10|1)(0|0)(0|00)
205211
int3c2e_envs = int3c2e_opt.int3c2e_envs
@@ -208,10 +214,12 @@ def proc(aux_batch_iter):
208214
buf = cp.empty((nao_pair*batch_size))
209215
buf2 = cp.empty((blksize, nao, nao))
210216
buf1 = cp.empty((blksize, nao, nocc))
211-
for kbatch in aux_batch_iter:
212-
if kbatch >= num_aux_batches:
217+
# Distribute packed auxiliary tiles through the existing device iterator.
218+
for batch_id in aux_batch_iter:
219+
if batch_id >= num_aux_batches:
213220
break
214-
naux_in_batch = aux_offsets[kbatch+1] - aux_offsets[kbatch]
221+
kbatch, kend = aux_splits[batch_id:batch_id+2]
222+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
215223
aux0 = aux1 = aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
216224
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf)
217225
for k0, k1 in lib.prange(0, naux_in_batch, blksize):
@@ -237,7 +245,7 @@ def proc(aux_batch_iter):
237245
ctypes.byref(int3c2e_envs),
238246
ctypes.c_int(shm_size_max),
239247
ctypes.c_int(len(shl_pair_offsets) - 1),
240-
ctypes.c_int(1),
248+
ctypes.c_int(kend-kbatch),
241249
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
242250
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
243251
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),

‎gpu4pyscf/df/hessian/uhf.py‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
from pyscf import lib
2828
from pyscf.gto import ATOM_OF
2929
from gpu4pyscf.lib import logger
30+
from gpu4pyscf.lib.utils import splits_by_blocksize
3031
from gpu4pyscf.lib.cupy_helper import (
3132
contract, asarray, ndarray, transpose_sum, get_avail_mem, fill_symmetric)
3233
from gpu4pyscf.df.int3c2e_bdiv import (
@@ -174,7 +175,6 @@ def proc():
174175
ksh_offsets_cpu = l_ctr_aux_offsets
175176
ksh_offsets_gpu = cp.asarray(ksh_offsets_cpu+mol.nbas, dtype=np.int32)
176177
aux_offsets = aux_loc[ksh_offsets_cpu]
177-
num_aux_batches = len(aux_offsets) - 1
178178

179179
if j_factor != 0:
180180
dm = contract('npi,nqi->pq', dm_factor_l, dm_factor_r)
@@ -187,8 +187,9 @@ def proc():
187187
buf = cp.empty((nao_pair*batch_size))
188188
buf2 = cp.empty((blksize, nao, nao))
189189
buf1 = cp.empty((blksize, nao, nocc))
190-
for kbatch in range(num_aux_batches):
191-
naux_in_batch = aux_offsets[kbatch+1] - aux_offsets[kbatch]
190+
aux_splits = splits_by_blocksize(aux_offsets, batch_size)
191+
for kbatch, kend in zip(aux_splits[:-1], aux_splits[1:]):
192+
naux_in_batch = aux_offsets[kend] - aux_offsets[kbatch]
192193
aux_ao_offset = aux_loc[ksh_offsets_cpu[kbatch]]
193194
compressed = ndarray((nao_pair, naux_in_batch), buffer=buf)
194195
for k0, k1 in lib.prange(0, naux_in_batch, blksize):
@@ -214,7 +215,7 @@ def proc():
214215
ctypes.byref(int3c2e_envs),
215216
ctypes.c_int(shm_size_max),
216217
ctypes.c_int(len(shl_pair_offsets) - 1),
217-
ctypes.c_int(1),
218+
ctypes.c_int(kend-kbatch),
218219
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
219220
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
220221
ctypes.cast(ksh_offsets_gpu[kbatch:].data.ptr, ctypes.c_void_p),

‎gpu4pyscf/df/int3c2e_bdiv.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -216,7 +216,7 @@ def int3c2e_evaluator(self, ao_pair_batch_size=None, aux_batch_size=None,
216216
l_ctr_aux_offsets, uniq_l_ctr_aux = _split_l_ctr_pattern(
217217
l_ctr_aux_offsets, uniq_l_ctr_aux, aux_batch_size)
218218
ksh_offsets_cpu = l_ctr_aux_offsets
219-
aux_splits = range(len(ksh_offsets_cpu))
219+
aux_splits = splits_by_blocksize(aux_loc[ksh_offsets_cpu], aux_batch_size)
220220
aux_offsets = aux_loc[ksh_offsets_cpu[aux_splits]]
221221
if reorder_aux:
222222
aux_sorting = argsort_aux(l_ctr_aux_offsets, uniq_l_ctr_aux)
@@ -800,7 +800,7 @@ def _int3c2e_ip1_evaluator(int3c2e_opt, scheme, batch_size,
800800

801801
ksh_offsets_cpu = l_ctr_aux_offsets
802802
ksh_offsets_gpu = cp.asarray(ksh_offsets_cpu+mol.nbas, dtype=np.int32)
803-
aux_splits = range(len(ksh_offsets_cpu))
803+
aux_splits = splits_by_blocksize(aux_loc[ksh_offsets_cpu], batch_size)
804804
aux_offsets = aux_loc[ksh_offsets_cpu[aux_splits]]
805805
kern = getattr(libvhf_rys, kern)
806806
int3c2e_envs = int3c2e_opt.int3c2e_envs
@@ -823,7 +823,7 @@ def evaluate_j3c(batch_id, out=None):
823823
ctypes.c_double(lr_factor), ctypes.c_double(sr_factor),
824824
ctypes.c_int(shm_size_max),
825825
ctypes.c_int(len(shl_pair_offsets) - 1),
826-
ctypes.c_int(1),
826+
ctypes.c_int(aux_split1 - aux_split0),
827827
ctypes.cast(shl_pair_offsets.data.ptr, ctypes.c_void_p),
828828
ctypes.cast(bas_ij_idx.data.ptr, ctypes.c_void_p),
829829
ctypes.cast(ksh_offsets_gpu[aux_split0:].data.ptr, ctypes.c_void_p),

0 commit comments

Comments
 (0)