Fix ACA workspace allocation for strain and max_iter=0 edge case

cutde/aca.cu

@@ -459,7 +459,7 @@ void aca_${name}(
     }
 
     if (team_idx == 0) {
-        n_terms[block_idx] = k + 1;
+        n_terms[block_idx] = max_iter > 0 ? k + 1 : 0;
     }
 }
 </%def>

cutde/aca.py

@@ -78,7 +78,7 @@ def call_clu_aca(
         gpu_buffer = backend.empty(block_sizes.sum(), float_type)
 
         # Storage for temporary rows and columns: RIref, RJref, RIstar, RJstar
-        fworkspace_per_block = n_cols + n_rows + 3 * n_cols + vec_dim * n_rows
+        fworkspace_per_block = n_cols + n_rows + vec_dim * n_cols + 3 * n_rows
         fworkspace_ends = np.cumsum(fworkspace_per_block)
         fworkspace_starts = fworkspace_ends - fworkspace_per_block
         gpu_fworkspace = backend.empty(fworkspace_ends[-1], float_type)

tests/test_aca.py

@@ -150,5 +150,73 @@ def runner(
         assert diff_frob < 1e-3
 
 
+@pytest.mark.parametrize("field_spec", [FS.DISP_SPEC, FS.STRAIN_SPEC])
+def test_aca_max_iter_zero(field_spec):
+    """When max_iter=0, no ACA terms should be produced and U @ V should be
+    the zero matrix."""
+    _, vec_dim = field_spec
+    n_obs, n_src = 5, 5
+    pts, tris, _ = setup_matrix_test(np.float64, False, n_obs=n_obs, n_src=n_src)
+
+    result = call_clu_aca(
+        pts,
+        tris,
+        [0],
+        [n_obs],
+        [0],
+        [n_src],
+        0.25,
+        [1e-4],
+        [0],
+        field_spec,
+    )
+
+    U, V = result[0]
+    assert U.shape == (n_obs * vec_dim, 0), f"Expected 0 terms, got U.shape={U.shape}"
+    assert V.shape == (0, n_src * 3), f"Expected 0 terms, got V.shape={V.shape}"
+    np.testing.assert_array_equal(U @ V, np.zeros((n_obs * vec_dim, n_src * 3)))
+
+
+def test_aca_strain_asymmetric_block():
+    """Regression test for fworkspace_per_block bug (gh-52).
+
+    With vec_dim=6 (strain) and n_src > 2*n_obs, the workspace was
+    undersized causing buffer overflow into adjacent blocks' memory."""
+    n_obs, n_src = 5, 20
+    pts, tris, _ = setup_matrix_test(np.float64, False, n_obs=n_obs, n_src=n_src)
+    pts[:, 0] -= 150
+
+    M = FS.strain_matrix(pts, tris, 0.25)
+    M = M.reshape((M.shape[0] * M.shape[1], M.shape[2] * M.shape[3]))
+
+    obs_starts = [0, 0]
+    obs_ends = [n_obs, n_obs]
+    src_starts = [0, 0]
+    src_ends = [n_src, n_src]
+
+    result = call_clu_aca(
+        pts,
+        tris,
+        obs_starts,
+        obs_ends,
+        src_starts,
+        src_ends,
+        0.25,
+        [1e-4] * 2,
+        [200] * 2,
+        FS.STRAIN_SPEC,
+        Iref0=np.zeros(2, dtype=np.int32),
+        Jref0=np.zeros(2, dtype=np.int32),
+    )
+
+    for i in range(2):
+        U, V = result[i]
+        diff = M - U @ V
+        diff_frob = np.sqrt(np.sum(diff**2))
+        assert (
+            diff_frob < 1e-3
+        ), f"Block {i}: Frobenius error {diff_frob:.3e} exceeds tolerance"
+
+
 if __name__ == "__main__":
     runner(np.float32, False, "disp", 40, compare_against_py=False, benchmark_iters=2)