Skip to content

Commit ffa2f72

Browse files
cheerly-pkucheerly-pku
committed
Keep PPCG block subspace two-block only
1 parent d31d9e2 commit ffa2f72

3 files changed

Lines changed: 17 additions & 86 deletions

File tree

‎source/source_hsolver/diago_ppcg.h‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -95,9 +95,6 @@ class DiagoPPCG
9595
std::vector<T> w_; // residual / preconditioned residual
9696
std::vector<T> sw_; // S * w
9797
std::vector<T> hw_; // H * w
98-
std::vector<T> p_; // previous search direction (for block subspace)
99-
std::vector<T> sp_; // S * p
100-
std::vector<T> hp_; // H * p
10198
std::vector<T> rr_psi_; // Rayleigh-Ritz rotation workspace
10299
std::vector<T> rr_spsi_;
103100
std::vector<T> rr_hpsi_;
@@ -165,7 +162,6 @@ class DiagoPPCG
165162
std::vector<T> m; // M matrix (projected S)
166163
std::vector<Real> eval; // eigenvalues
167164
std::vector<Real> w_scale;
168-
std::vector<Real> p_scale;
169165
};
170166

171167
void lock_epairs(const std::vector<T>& residual,
@@ -174,15 +170,13 @@ class DiagoPPCG
174170

175171
void build_small_subspace(const T* psi,
176172
const std::vector<int>& cols,
177-
bool use_p,
178173
SmallSubspace& subspace) const;
179174

180175
void solve_small_generalized(int dim, SmallSubspace& subspace) const;
181176

182177
void update_one_block(T* psi,
183178
const std::vector<int>& cols,
184179
int l,
185-
bool use_p,
186180
const SmallSubspace& subspace);
187181

188182
bool is_s_orthonormal(const T* psi, const T* spsi, int ncol) const;

‎source/source_hsolver/ppcg/diago_ppcg_diag.hpp‎

Lines changed: 7 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,6 @@ double DiagoPPCG<T, Device>::diag(const HPsiFunc& hpsi_func,
3030
w_.assign(sz, T(0));
3131
sw_.assign(sz, T(0));
3232
hw_.assign(sz, T(0));
33-
p_.clear();
34-
sp_.clear();
35-
hp_.clear();
3633
rr_psi_.resize(sz);
3734
rr_spsi_.resize(sz);
3835
rr_hpsi_.resize(sz);
@@ -118,15 +115,10 @@ double DiagoPPCG<T, Device>::diag(const HPsiFunc& hpsi_func,
118115

119116
avg_iter += static_cast<double>(nact) / static_cast<double>(ncol);
120117

121-
// Use the stable 2-block [psi, w] projected subspace.
122-
// The historical p block is kept in the implementation helpers,
123-
// but is not enabled in the production path because it can make
124-
// the small generalized eigenproblem indefinite on common test
125-
// cases.
126-
// w is normalized to unit S-norm before building the
127-
// Gram matrix (see build_small_subspace), which keeps M
128-
// well-conditioned even when residuals are small.
129-
const bool use_p_now = false;
118+
// Use the stable 2-block [psi, w] projected subspace. The
119+
// preconditioned residual w is normalized to unit S-norm before
120+
// building the Gram matrix (see build_small_subspace), which
121+
// keeps M well-conditioned even when residuals are small.
130122

131123
// Block subspace solve.
132124
for (int isb = 0; isb < nsb; ++isb)
@@ -136,9 +128,9 @@ double DiagoPPCG<T, Device>::diag(const HPsiFunc& hpsi_func,
136128
cols.assign(active_cols.begin() + i0,
137129
active_cols.begin() + i0 + l);
138130

139-
build_small_subspace(psi_in, cols, use_p_now, subspace);
140-
solve_small_generalized((use_p_now ? 3 : 2) * l, subspace);
141-
update_one_block(psi_in, cols, l, use_p_now, subspace);
131+
build_small_subspace(psi_in, cols, subspace);
132+
solve_small_generalized(2 * l, subspace);
133+
update_one_block(psi_in, cols, l, subspace);
142134
}
143135

144136
// Rayleigh-Ritz after each block update keeps the global subspace

‎source/source_hsolver/ppcg/diago_ppcg_subspace.hpp‎

Lines changed: 10 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -37,46 +37,36 @@ void DiagoPPCG<T, Device>::lock_epairs(
3737
}
3838

3939
// ---------------------------------------------------------------------------
40-
// Build K = V^H H V and M = V^H S V where V = [psi, w, p]
40+
// Build K = V^H H V and M = V^H S V where V = [psi, w]
4141
// ---------------------------------------------------------------------------
4242
template <typename T, typename Device>
4343
void DiagoPPCG<T, Device>::build_small_subspace(
4444
const T* psi,
4545
const std::vector<int>& cols,
46-
bool use_p,
4746
SmallSubspace& subspace) const
4847
{
4948
const int l = static_cast<int>(cols.size());
50-
const int nblk = use_p ? 3 : 2;
51-
const int dim = nblk * l;
49+
const int dim = 2 * l;
5250
subspace.k.resize(dim * dim);
5351
subspace.m.resize(dim * dim);
5452
subspace.eval.resize(dim);
5553
subspace.w_scale.assign(l, static_cast<Real>(1));
56-
subspace.p_scale.assign(l, static_cast<Real>(1));
5754

5855
std::vector<T> psi_l, spsi_l, hpsi_l;
5956
std::vector<T> w_l, sw_l, hw_l;
60-
std::vector<T> p_l, sp_l, hp_l;
6157
copy_cols(psi, cols, psi_l);
6258
copy_cols(spsi_.data(), cols, spsi_l);
6359
copy_cols(hpsi_.data(), cols, hpsi_l);
6460
copy_cols(w_.data(), cols, w_l);
6561
copy_cols(sw_.data(), cols, sw_l);
6662
copy_cols(hw_.data(), cols, hw_l);
67-
if (use_p)
68-
{
69-
copy_cols(p_.data(), cols, p_l);
70-
copy_cols(sp_.data(), cols, sp_l);
71-
copy_cols(hp_.data(), cols, hp_l);
72-
}
7363

7464
// ---------------------------------------------------------------------------
75-
// Normalize w and p columns to unit S-norm for numerical stability.
65+
// Normalize w columns to unit S-norm for numerical stability.
7666
//
77-
// The [w, p] block of the Gram matrix M has entries O(||w||²) which
78-
// become tiny when residuals are small, making M nearly singular and
79-
// causing sygvd to produce garbage eigenvectors.
67+
// The w block of the Gram matrix M has entries O(||w||^2) which become
68+
// tiny when residuals are small, making M nearly singular and causing
69+
// sygvd to produce garbage eigenvectors.
8070
//
8171
// Scaling to unit S-norm keeps M well-conditioned (diagonal ~1) without
8272
// changing the subspace. The Ritz values are identical and the Ritz
@@ -117,8 +107,6 @@ void DiagoPPCG<T, Device>::build_small_subspace(
117107
}
118108
};
119109
scale_to_unit_snorm(w_l, sw_l, hw_l, l, subspace.w_scale);
120-
if (use_p)
121-
scale_to_unit_snorm(p_l, sp_l, hp_l, l, subspace.p_scale);
122110

123111
auto copy_block = [&](const std::vector<T>& src,
124112
const int col0,
@@ -157,12 +145,6 @@ void DiagoPPCG<T, Device>::build_small_subspace(
157145
copy_block(w_l, l, basis);
158146
copy_block(hw_l, l, hbasis);
159147
copy_block(sw_l, l, sbasis);
160-
if (use_p)
161-
{
162-
copy_block(p_l, 2 * l, basis);
163-
copy_block(hp_l, 2 * l, hbasis);
164-
copy_block(sp_l, 2 * l, sbasis);
165-
}
166148

167149
gram(basis.data(), hbasis.data(), dim, dim, subspace.k, dim);
168150
gram(basis.data(), sbasis.data(), dim, dim, subspace.m, dim);
@@ -225,36 +207,25 @@ void DiagoPPCG<T, Device>::update_one_block(
225207
T* psi,
226208
const std::vector<int>& cols,
227209
int l,
228-
bool use_p,
229210
const SmallSubspace& subspace)
230211
{
231-
const int dim = (use_p ? 3 : 2) * l;
212+
const int dim = 2 * l;
232213
const T* eigvec = subspace.k.data();
233214

234215
std::vector<T> psi_l, spsi_l, hpsi_l;
235216
std::vector<T> w_l, sw_l, hw_l;
236-
std::vector<T> p_l, sp_l, hp_l;
237217
copy_cols(psi, cols, psi_l);
238218
copy_cols(spsi_.data(), cols, spsi_l);
239219
copy_cols(hpsi_.data(), cols, hpsi_l);
240220
copy_cols(w_.data(), cols, w_l);
241221
copy_cols(sw_.data(), cols, sw_l);
242222
copy_cols(hw_.data(), cols, hw_l);
243-
if (use_p)
244-
{
245-
copy_cols(p_.data(), cols, p_l);
246-
copy_cols(sp_.data(), cols, sp_l);
247-
copy_cols(hp_.data(), cols, hp_l);
248-
}
249223

250224
std::vector<T> psi_new(ld_psi_ * l, T(0));
251225
std::vector<T> spsi_new(ld_psi_ * l, T(0));
252226
std::vector<T> hpsi_new(ld_psi_ * l, T(0));
253227

254228
std::vector<T> coeff_state(dim * l, T(0));
255-
std::vector<T> coeff_dir;
256-
if (use_p)
257-
coeff_dir.assign(dim * l, T(0));
258229
#ifdef _OPENMP
259230
#pragma omp parallel for schedule(static) if (l * l > 4096)
260231
#endif
@@ -265,19 +236,11 @@ void DiagoPPCG<T, Device>::update_one_block(
265236
coeff_state[i + j * dim] = eigvec[i + j * dim];
266237
const T cw = eigvec[(l + i) + j * dim] * subspace.w_scale[i];
267238
coeff_state[(l + i) + j * dim] = cw;
268-
if (use_p)
269-
{
270-
coeff_dir[(l + i) + j * dim] = cw;
271-
const T cp = eigvec[(2*l + i) + j * dim] * subspace.p_scale[i];
272-
coeff_state[(2*l + i) + j * dim] = cp;
273-
coeff_dir[(2*l + i) + j * dim] = cp;
274-
}
275239
}
276240
}
277241

278242
auto fill_basis = [&](const std::vector<T>& a,
279243
const std::vector<T>& b,
280-
const std::vector<T>& c,
281244
std::vector<T>& basis)
282245
{
283246
basis.resize(ld_psi_ * dim);
@@ -292,12 +255,6 @@ void DiagoPPCG<T, Device>::update_one_block(
292255
std::copy(b.begin() + j * ld_psi_,
293256
b.begin() + (j + 1) * ld_psi_,
294257
basis.begin() + (l + j) * ld_psi_);
295-
if (use_p)
296-
{
297-
std::copy(c.begin() + j * ld_psi_,
298-
c.begin() + (j + 1) * ld_psi_,
299-
basis.begin() + (2 * l + j) * ld_psi_);
300-
}
301258
}
302259
};
303260

@@ -325,9 +282,9 @@ void DiagoPPCG<T, Device>::update_one_block(
325282
std::vector<T> psi_basis;
326283
std::vector<T> spsi_basis;
327284
std::vector<T> hpsi_basis;
328-
fill_basis(psi_l, w_l, p_l, psi_basis);
329-
fill_basis(spsi_l, sw_l, sp_l, spsi_basis);
330-
fill_basis(hpsi_l, hw_l, hp_l, hpsi_basis);
285+
fill_basis(psi_l, w_l, psi_basis);
286+
fill_basis(spsi_l, sw_l, spsi_basis);
287+
fill_basis(hpsi_l, hw_l, hpsi_basis);
331288

332289
combine(psi_basis, coeff_state, psi_new);
333290
combine(spsi_basis, coeff_state, spsi_new);
@@ -336,18 +293,6 @@ void DiagoPPCG<T, Device>::update_one_block(
336293
scatter_cols(psi, cols, psi_new);
337294
scatter_cols(spsi_.data(), cols, spsi_new);
338295
scatter_cols(hpsi_.data(), cols, hpsi_new);
339-
if (use_p)
340-
{
341-
std::vector<T> p_new(ld_psi_ * l, T(0));
342-
std::vector<T> sp_new(ld_psi_ * l, T(0));
343-
std::vector<T> hp_new(ld_psi_ * l, T(0));
344-
combine(psi_basis, coeff_dir, p_new);
345-
combine(spsi_basis, coeff_dir, sp_new);
346-
combine(hpsi_basis, coeff_dir, hp_new);
347-
scatter_cols(p_.data(), cols, p_new);
348-
scatter_cols(sp_.data(), cols, sp_new);
349-
scatter_cols(hp_.data(), cols, hp_new);
350-
}
351296
}
352297

353298
} // namespace hsolver

0 commit comments

Comments
 (0)