[DiadkoShmek]

TriangularLinearOperator.solve forwards left_tensor into delegated solves for repeated or non-dense representations, then multiplies by it again. A square left factor produces L @ L @ solve(A, R) instead of L @ solve(A, R); a rectangular factor can raise a shape error. This also affects the quadratic-form case with L = R.mT.

Let the delegated solve compute only solve(A, R). Apply the wrapper's left factor once while a vector RHS still has its normalized one-column shape, then restore vector output shape.

For a public triangular-plus-diagonal construction with A = diag(4, 6), R = [[1], [3]], and L = [[2, 1], [0, 3]], the current result is [[3.5], [4.5]]; the corrected and exact result is [[1], [1.5]].

Validation on PyTorch 2.14.1 CPU, float32/float64:

• The three added unittest methods cover 36 subcases. On the base these produce 12 failures and 24 errors; the candidate passes all 127 selected lower/upper triangular unittest methods, including the new subcases.
• A separate NumPy analytic value/first-gradient oracle passes 100 configurations. Independent review passes 192 checks (120 base failures), including 24 no-left controls and 12 tied-RHS quadratic-form cases.
• The changed files pass the configured flake8 5.0.4 and ufmt 2.0.0 checks. Source imports and actual Torch library hashes were bound before and after execution.

The checks cover added-diagonal SPD examples and supported repeated general triangles, shared RHS/vector forms and selected broadcasting. Additional RHS batch prefixes retain a paired pre-existing BatchRepeat limitation. The full repository suite, GPU and higher-order AD were not tested.

I directed this work and take responsibility for the contribution and its follow-up. AI assistance was used in producing the change.