Add damp support to CG solver (closes #406)

[gaoflow]

Motivation

CG currently only solves $\mathbf{Op},\mathbf{x} = \mathbf{y}$. As discussed in #406, it should also be able to solve the damped system

$$(\mathbf{Op} + \epsilon\mathbf{I}),\mathbf{x} = \mathbf{y},$$

consistent with cgls/lsqr (which both expose damp). The damp argument used to exist on CG as an unused stub and was removed in v2; this PR re-introduces it with a working implementation.

This matches the issue's definition of done ("Added implementation with damp != 0") and the semantics @mrava87 confirmed there: solve $(A + \text{damp},I)x = y$, equivalent to adding damp*I to the operator and using damp=0.

Implementation

Damping is applied by adding $\epsilon\mathbf{c}$ to every application of the operator $\mathbf{Op},\mathbf{c}$ in the iterations (and $\epsilon\mathbf{x}_0$ to the initial residual when x0 is given) — i.e. the operator is implicitly replaced by $\mathbf{Op} + \epsilon\mathbf{I}$:

Opc = self.Op.matvec(self.c)
if self.damp != 0.0:
    Opc = Opc + self.damp * self.c

damp is threaded through cg(), CG.solve and CG.setup (default 0.0). With damp=0.0 the code path is identical to before, so existing behaviour is unchanged. Op must remain square and SPD, which $\mathbf{Op} + \epsilon\mathbf{I}$ is for $\epsilon \ge 0$.

Verification

Added test_cg_damp which checks, for both preallocate modes, the equivalence @mrava87 suggested in the issue:

• internal damping cg(Op, y, damp=d) == adding d*I to the operator cg(Op + d*I, y, damp=0)
• both == dense (Op + d*I)^{-1} y
• damp=0 still recovers the undamped Op^{-1} y

Full pytests/test_solver.py passes (100 tests), ruff check clean.

[codacy-production]

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 complexity · 6 duplication

Metric	Results
Complexity	0
Duplication	6

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer
TIP This summary will be updated as you push new changes.}

[mrava87]

@gaoflow thanks for a nice contribution :)

I left some minor comments, mostly to ensure consistency with other bits of codes in other solvers having the damp option.

[gaoflow]

Thanks for the review! Addressed all points in 6828c57:

• damp docstrings shortened to Damping coefficient in CG.setup, CG.solve and cg, and the CG Notes reduced to a single damped-system sentence.
• Initial damping correction in setup now follows the non-preallocate (self.r = self.r - damp * x) vs preallocate (ncp.subtract(..., out=self.r)) pattern used elsewhere.
• Trimmed the step comment to # add damping contribution.
• test_cg_damp: now builds y = Aop * x from a known x, renamed the loop variable to preallocate, and dropped the issue ref from the docstring.

(The red build check is the self-hosted GPU runner failing on a ptxas version mismatch — Unsupported .version 8.8 — in unrelated test_kirchhoff/test_nonstatconvolve CUDA tests; all matrix, Lint and Azure builds are green.)

[mrava87]

Thanks for the review! Addressed all points in 6828c57:

• damp docstrings shortened to Damping coefficient in CG.setup, CG.solve and cg, and the CG Notes reduced to a single damped-system sentence.
• Initial damping correction in setup now follows the non-preallocate (self.r = self.r - damp * x) vs preallocate (ncp.subtract(..., out=self.r)) pattern used elsewhere.
• Trimmed the step comment to # add damping contribution.
• test_cg_damp: now builds y = Aop * x from a known x, renamed the loop variable to preallocate, and dropped the issue ref from the docstring.

(The red build check is the self-hosted GPU runner failing on a ptxas version mismatch — Unsupported .version 8.8 — in unrelated test_kirchhoff/test_nonstatconvolve CUDA tests; all matrix, Lint and Azure builds are green.)

Thanks. Please reply to each comment with link to commit you fixed the issue.