CLAUDE: rung C1 - resident prepareLMA (assemble + serial normalize on-device)

Design 3-A4i rung C1: the LMA measured inputs never leave the GPU.
- pose_lma_assemble (grid over resident corr rows): peaks {dx,dy,str,e0x,
  e0y,l0,l1,valid} + packed corr indices -> per-tile y offsets, raw strength
  weights, eigen 2x2 transform (setEigenTransform math, float; identity when
  use_eigen==0) and the fx selection mask. Parameterized contract (counts,
  slot, shift, components, policy scalars) - CUAS specifics stay host policy.
- pose_lma_prepare_norm + pose_lma_prepare_finish: SERIAL one-thread tails
  (sample normalize + pull rows; y+=fx, H diagonal in Java getWJtJlambda
  order, regularization weights, full-weight normalize + pure_weight) -
  serial by design so the Java/C++ float clones gate BIT-EXACT.
- tp_proc_exec_pose_lma_prepare_resident: camera/pose/centers/policy H2D,
  assemble -> norm -> pose_fx_jacobian(current, NO readback) -> finish;
  D2H = 8-float compact result (+optional capture of y/weights/eigen for
  the one-shot oracle). Same scratch bookkeeping as the step so a following
  no-H2D step never regrows.
- tp_proc_exec_pose_lma_resident_step: all-NULL weights/y/eigen = use the
  resident prepared buffers (skip H2D); mixed NULLs rejected.
Gates: direct CUDA prepare 5120/159 BIT-EXACT vs serial float oracle
(y/weights/eigen/selection/scalars all 0 mismatches); resident + 19-float +
fx_jacobian NVRTC PASS; pose_corr bit-exact @tol 0; compute-sanitizer 0
errors on both binaries.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

src/tests/test_pose_lma_step.cu

@@ -293,5 +293,145 @@ int main(int argc, char ** argv)
 			preparedBad, maxProducts, maxDelta, maxCandidate, maxRms,
 			residentOk ? "PASS" : "FAIL");
 	if (!residentOk) bad++;
+
+	// ── Rung C1: resident prepare kernels vs serial float CPU oracle ─────────
+	// assemble is a per-tile scatter (order-free); both tails are SERIAL on the
+	// GPU by design, so the whole prepare is BIT-EXACT against this CPU clone.
+	// By Claude on 07/17/2026.
+	{
+		const int numTilesP = 5120, numCorr = 200, slot = 0xff, shift = 8;
+		const int nsP = POSE_NUM_COMPONENTS * numTilesP;
+		const int numValuesP = nsP + POSE_NUM_PARAMS;
+		const float eigMin = 0.5f, eigMax = 20.0f, minConf = 0.0f;
+		std::vector<int> corrIdx(numCorr);
+		std::vector<float> peaks((size_t) numCorr * 8, 0.0f);
+		unsigned int lcg = 12345u;
+		auto nextf = [&lcg](float lo, float hi) {
+			lcg = lcg * 1664525u + 1013904223u;
+			return lo + (hi - lo) * (float)(lcg >> 8) / 16777216.0f; };
+		std::vector<float> rawFxP((size_t) nsP, 0.0f);
+		std::vector<float> rawJtP((size_t) POSE_NUM_PARAMS * nsP, 0.0f);
+		std::vector<unsigned char> validP(numTilesP, 0);
+		for (int n = 0; n < numCorr; n++) {
+			const int tile = 30 + 25 * n; // distinct tiles, < 5120
+			const bool otherSlot = (n % 17) == 3;      // wrong pair slot -> skipped
+			corrIdx[n] = (tile << shift) | (otherSlot ? 0x03 : slot);
+			float* row = &peaks[(size_t) n * 8];
+			row[0] = nextf(-2.0f, 2.0f);               // dx
+			row[1] = nextf(-2.0f, 2.0f);               // dy
+			row[2] = ((n % 13) == 5) ? -0.5f : nextf(0.2f, 2.0f); // some non-positive str
+			const float ang = nextf(0.0f, 6.28f);
+			row[3] = std::cos(ang); row[4] = std::sin(ang);       // e0
+			row[5] = nextf(0.3f, 30.0f); row[6] = nextf(0.3f, 30.0f); // l0,l1 (some < eigMin^2)
+			row[7] = ((n % 11) == 7) ? 0.0f : 1.0f;    // some kernel-rejected
+			const int vi = POSE_NUM_COMPONENTS * tile;
+			validP[tile] = ((n % 19) == 9) ? 0 : 1;    // some invalid raw projections
+			rawFxP[vi] = nextf(-1.0f, 1.0f); rawFxP[vi + 1] = nextf(-1.0f, 1.0f);
+			for (int par = 0; par < POSE_NUM_PARAMS; par++) {
+				rawJtP[(size_t) par * nsP + vi] = nextf(-300.0f, 300.0f);
+				rawJtP[(size_t) par * nsP + vi + 1] = nextf(-300.0f, 300.0f);
+			}
+		}
+		const float pullT[3] = {0.01f, -0.02f, 0.003f};
+		const float regW[3] = {0.05f, 0.05f, 0.02f};
+		// CPU serial float oracle (mirrors IntersceneLmaFloat.prepareResidentOracle)
+		std::vector<float> expY(numValuesP, 0.0f), expW(numValuesP, 0.0f);
+		std::vector<float> expE((size_t) 4 * numTilesP, 0.0f);
+		std::vector<unsigned char> expSel(numTilesP, 0);
+		for (int n = 0; n < numCorr; n++) {
+			if ((corrIdx[n] & ((1 << shift) - 1)) != slot) continue;
+			const int tile = corrIdx[n] >> shift;
+			const float* row = &peaks[(size_t) n * 8];
+			if (row[7] == 0.0f) continue;
+			const float str = row[2];
+			if (!(str >= minConf) || !(str > 0.0f)) continue;
+			const int vi = POSE_NUM_COMPONENTS * tile, ei = 4 * tile;
+			expSel[tile] = 1; expW[vi] = str;
+			expY[vi] = row[0]; expY[vi + 1] = row[1];
+			const float k0inv = 1.0f / eigMax;
+			const float k0 = std::fmax(0.0f, 1.0f / std::fmax(eigMin, std::sqrt(row[5])) - k0inv);
+			const float k1 = std::fmax(0.0f, 1.0f / std::fmax(eigMin, std::sqrt(row[6])) - k0inv);
+			expE[ei] = k0 * row[3]; expE[ei + 1] = k0 * row[4];
+			expE[ei + 2] = k1 * row[4]; expE[ei + 3] = k1 * -row[3];
+		}
+		float sumRaw = 0.0f; int nMeas = 0;
+		for (int t = 0; t < numTilesP; t++) {
+			const float w = expW[POSE_NUM_COMPONENTS * t];
+			sumRaw += w + w; if (w > 0.0f) nMeas++; }
+		const float sN = 1.0f / sumRaw;
+		for (int t = 0; t < numTilesP; t++) { const int vi = POSE_NUM_COMPONENTS * t;
+			const float w = expW[vi] * sN; expW[vi] = w; expW[vi + 1] = w; }
+		for (int par = 0; par < POSE_NUM_PARAMS; par++) expY[nsP + par] = pullT[par];
+		for (int t = 0; t < numTilesP; t++) { const int vi = POSE_NUM_COMPONENTS * t;
+			if ((expW[vi] > 0.0f) && validP[t]) {
+				expY[vi] = expY[vi] + rawFxP[vi]; expY[vi + 1] = expY[vi + 1] + rawFxP[vi + 1]; } }
+		float hD[3] = {0.0f, 0.0f, 0.0f};
+		for (int t = 0; t < numTilesP; t++) { const int vi = POSE_NUM_COMPONENTS * t;
+			const float w = expW[vi]; if (!(w > 0.0f)) continue;
+			const bool have = validP[t] != 0; const int ei = 4 * t;
+			for (int par = 0; par < POSE_NUM_PARAMS; par++) {
+				const float jx = have ? rawJtP[(size_t) par * nsP + vi] : 0.0f;
+				const float jy = have ? rawJtP[(size_t) par * nsP + vi + 1] : 0.0f;
+				const float tj0 = expE[ei] * jx + expE[ei + 1] * jy;
+				const float tj1 = expE[ei + 2] * jx + expE[ei + 3] * jy;
+				hD[par] = hD[par] + w * (tj0 * tj0);
+				hD[par] = hD[par] + w * (tj1 * tj1); } }
+		float expReg[3];
+		for (int par = 0; par < POSE_NUM_PARAMS; par++) {
+			expReg[par] = (hD[par] > 0.0f) ? (regW[par] / std::sqrt(hD[par])) : 0.0f;
+			expW[nsP + par] = expReg[par]; }
+		float pureS = 0.0f; for (int i = 0; i < nsP; i++) pureS += expW[i];
+		float fullS = pureS; for (int par = 0; par < POSE_NUM_PARAMS; par++) fullS += expW[nsP + par];
+		const float scaleS = 1.0f / fullS;
+		for (int i = 0; i < numValuesP; i++) expW[i] *= scaleS;
+		// GPU
+		GpuBuf<int> gIdx(numCorr); GpuBuf<float> gPeaks(peaks.size());
+		GpuBuf<float> gY(numValuesP), gW(numValuesP), gE((size_t) 4 * numTilesP);
+		GpuBuf<unsigned char> gSel(numTilesP), gValid(numTilesP);
+		GpuBuf<float> gFx(rawFxP.size()), gJt(rawJtP.size());
+		GpuBuf<float> gPull(3), gReg(3), gPrep(POSE_LMA_PREP_RESULT);
+		gIdx.upload(corrIdx.data()); gPeaks.upload(peaks.data());
+		gValid.upload(validP.data()); gFx.upload(rawFxP.data()); gJt.upload(rawJtP.data());
+		gPull.upload(pullT); gReg.upload(regW);
+		checkCudaErrors(cudaMemset(gY.dev(), 0, numValuesP * sizeof(float)));
+		checkCudaErrors(cudaMemset(gW.dev(), 0, numValuesP * sizeof(float)));
+		checkCudaErrors(cudaMemset(gE.dev(), 0, (size_t) 4 * numTilesP * sizeof(float)));
+		checkCudaErrors(cudaMemset(gSel.dev(), 0, numTilesP));
+		checkCudaErrors(cudaMemset(gPrep.dev(), 0, POSE_LMA_PREP_RESULT * sizeof(float)));
+		pose_lma_assemble<<<(numCorr + 255) / 256, 256>>>(numCorr, numTilesP, slot, shift,
+				POSE_NUM_COMPONENTS, 1, eigMin, eigMax, minConf, 0,
+				gIdx.dev(), gPeaks.dev(), gY.dev(), gW.dev(), gE.dev(), gSel.dev());
+		pose_lma_prepare_norm<<<1, 1>>>(numTilesP, POSE_NUM_COMPONENTS, POSE_NUM_PARAMS,
+				gPull.dev(), gY.dev(), gW.dev(), gPrep.dev());
+		pose_lma_prepare_finish<<<1, 1>>>(numTilesP, POSE_NUM_COMPONENTS, POSE_NUM_PARAMS,
+				gFx.dev(), gJt.dev(), gValid.dev(), gE.dev(), gReg.dev(),
+				gY.dev(), gW.dev(), gPrep.dev());
+		checkCudaErrors(cudaGetLastError()); checkCudaErrors(cudaDeviceSynchronize());
+		const std::vector<float> gotY = gY.download();
+		const std::vector<float> gotW = gW.download();
+		const std::vector<float> gotE = gE.download();
+		const std::vector<unsigned char> gotSel = gSel.download();
+		const std::vector<float> gotPrep = gPrep.download();
+		int badBits = 0, badSel = 0;
+		for (int i = 0; i < numValuesP; i++) {
+			if (std::memcmp(&gotY[i], &expY[i], sizeof(float)) != 0) badBits++;
+			if (std::memcmp(&gotW[i], &expW[i], sizeof(float)) != 0) badBits++; }
+		for (size_t i = 0; i < expE.size(); i++)
+			if (std::memcmp(&gotE[i], &expE[i], sizeof(float)) != 0) badBits++;
+		for (int t = 0; t < numTilesP; t++) if (gotSel[t] != expSel[t]) badSel++;
+		const bool prepOk = (badBits == 0) && (badSel == 0) &&
+				(gotPrep[POSE_LMA_PREP_VALID] == 1.0f) &&
+				(gotPrep[POSE_LMA_PREP_NUM_MEAS] == (float) nMeas) &&
+				(std::memcmp(&gotPrep[POSE_LMA_PREP_SUM_W_RAW], &sumRaw, sizeof(float)) == 0) &&
+				(std::memcmp(&gotPrep[POSE_LMA_PREP_FULL_WEIGHT], &fullS, sizeof(float)) == 0) &&
+				(std::memcmp(&gotPrep[POSE_LMA_PREP_REG_W0], &expReg[0], sizeof(float)) == 0) &&
+				(std::memcmp(&gotPrep[POSE_LMA_PREP_REG_W0 + 1], &expReg[1], sizeof(float)) == 0) &&
+				(std::memcmp(&gotPrep[POSE_LMA_PREP_REG_W0 + 2], &expReg[2], sizeof(float)) == 0);
+		std::printf("pose_lma prepare 5120/%d meas: bit mismatches=%d, selection mismatches=%d, "
+				"num_meas=%g/%d, pure_weight=%g -> %s\n",
+				nMeas, badBits, badSel, gotPrep[POSE_LMA_PREP_NUM_MEAS], nMeas,
+				gotPrep[POSE_LMA_PREP_PURE_WEIGHT], prepOk ? "PASS" : "FAIL");
+		if (!prepOk) bad++;
+	}
 	return bad ? EXIT_FAILURE : EXIT_SUCCESS;
 }

src/tp_lma.cu

@@ -769,6 +769,177 @@ extern "C" __global__ void pose_lma_reduce_rms(
 	}
 }
 
+// ── Rung C1: resident prepareLMA (design 2026-07-16) ─────────────────────────
+// See tp_lma.h for the contracts. The two tails are deliberately SERIAL
+// one-thread kernels: their loop order replicates the Java float-clone oracle
+// exactly (bit-exact gate), and O(num_values) float adds are microseconds.
+// By Claude on 07/17/2026.
+
+extern "C" __global__ void pose_lma_assemble(
+		int                   num_corr,
+		int                   num_tiles,
+		int                   corr_slot,
+		int                   ntile_shift,
+		int                   num_components,
+		int                   use_eigen,
+		float                 eig_min_sqrt,
+		float                 eig_max_sqrt,
+		float                 min_confidence,
+		int                   same_weights,
+		const int           * corr_indices,
+		const float         * peaks,
+		float               * y,
+		float               * weights,
+		float               * eigen,
+		unsigned char       * selection)
+{
+	const int n = (int) (blockIdx.x * blockDim.x + threadIdx.x);
+	if (n >= num_corr) return;
+	const int packed = corr_indices[n];
+	if ((packed & ((1 << ntile_shift) - 1)) != corr_slot) return;
+	const int tile = packed >> ntile_shift;
+	if ((tile < 0) || (tile >= num_tiles)) return;
+	const int row = n * 8; // PEAK_ROW_FLOATS {dx,dy,str,e0x,e0y,l0,l1,valid}
+	if (peaks[row + 7] == 0.0f) return;               // kernel-rejected peak
+	const float str = peaks[row + 2];
+	if (!(str >= min_confidence)) return;             // NaN strength fails here too
+	if (use_eigen) {
+		// Java setEigenTransform: a NaN eigenvector leaves transform null and
+		// setSamplesWeights then SKIPS the tile entirely. Non-finite lambdas are
+		// also skipped here (Java would poison the row with NaN; peaks from
+		// corr2D_peak_eig_f are finite for valid rows).
+		if (!isfinite(peaks[row + 3]) || !isfinite(peaks[row + 4]) ||
+		    !isfinite(peaks[row + 5]) || !isfinite(peaks[row + 6])) return;
+	}
+	if (!(str > 0.0f)) return;
+	const float w = same_weights ? 1.0f : str;
+	const int vi = num_components * tile;
+	selection[tile] = 1;
+	weights[vi] = w;                                   // raw; normalized by the tails
+	y[vi]     = peaks[row];                            // measured offsets; fx added
+	y[vi + 1] = peaks[row + 1];                        // by pose_lma_prepare_finish
+	const int ei = 4 * tile;
+	if (use_eigen) {
+		// setEigenTransform (float): am = {{e0x,e0y},{e0y,-e0x}}, row i scaled by
+		// k_i = max(0, 1/max(eig_min_sqrt, sqrt(l_i)) - 1/eig_max_sqrt).
+		const float e0x = peaks[row + 3];
+		const float e0y = peaks[row + 4];
+		const float k0inv = __fdiv_rn(1.0f, eig_max_sqrt);
+		const float k0 = fmaxf(0.0f, pose_lma_sub(
+				__fdiv_rn(1.0f, fmaxf(eig_min_sqrt, sqrtf(peaks[row + 5]))), k0inv));
+		const float k1 = fmaxf(0.0f, pose_lma_sub(
+				__fdiv_rn(1.0f, fmaxf(eig_min_sqrt, sqrtf(peaks[row + 6]))), k0inv));
+		eigen[ei]     = pose_lma_mul(k0, e0x);
+		eigen[ei + 1] = pose_lma_mul(k0, e0y);
+		eigen[ei + 2] = pose_lma_mul(k1, e0y);
+		eigen[ei + 3] = pose_lma_mul(k1, -e0x);
+	} else { // identity (Java flattenEigenTransform with null source)
+		eigen[ei]     = 1.0f;
+		eigen[ei + 3] = 1.0f;
+	}
+}
+
+extern "C" __global__ void pose_lma_prepare_norm(
+		int                   num_tiles,
+		int                   num_components,
+		int                   num_params,
+		const float         * pull_targets,
+		float               * y,
+		float               * weights,
+		float               * prep_result)
+{
+	if ((blockIdx.x != 0) || (threadIdx.x != 0)) return;
+	// Java setSamplesWeights: sum 2*w over measured tiles (w==0 adds exactly 0),
+	// scale by 1/sum, duplicate into the second component.
+	float sum = 0.0f;
+	int num_meas = 0;
+	for (int tile = 0; tile < num_tiles; tile++) {
+		const float w = weights[num_components * tile];
+		sum = pose_lma_add(sum, pose_lma_add(w, w));
+		if (w > 0.0f) num_meas++;
+	}
+	const float s = __fdiv_rn(1.0f, sum);
+	for (int tile = 0; tile < num_tiles; tile++) {
+		const int vi = num_components * tile;
+		const float w = pose_lma_mul(weights[vi], s);
+		weights[vi] = w;
+		weights[vi + 1] = w;
+	}
+	const int ns = num_components * num_tiles;
+	for (int par = 0; par < num_params; par++) {
+		y[ns + par] = pull_targets[par]; // pull-target rows (full RMS regularization)
+	}
+	prep_result[POSE_LMA_PREP_NUM_MEAS] = (float) num_meas;
+	prep_result[POSE_LMA_PREP_SUM_W_RAW] = sum;
+}
+
+extern "C" __global__ void pose_lma_prepare_finish(
+		int                   num_tiles,
+		int                   num_components,
+		int                   num_params,
+		const float         * raw_fx,
+		const float         * raw_jt,
+		const unsigned char * valid,
+		const float         * eigen,
+		const float         * param_regweights,
+		float               * y,
+		float               * weights,
+		float               * prep_result)
+{
+	if ((blockIdx.x != 0) || (threadIdx.x != 0)) return;
+	const int ns = num_components * num_tiles;
+	const int num_values = ns + num_params;
+	// y += fx(current) on measured tiles with valid raw projection (Java:
+	// y = fx.clone() + measured offsets; provider leaves invalid rows at 0).
+	for (int tile = 0; tile < num_tiles; tile++) {
+		const int vi = num_components * tile;
+		if ((weights[vi] > 0.0f) && valid[tile]) {
+			y[vi]     = pose_lma_add(y[vi],     raw_fx[vi]);
+			y[vi + 1] = pose_lma_add(y[vi + 1], raw_fx[vi + 1]);
+		}
+	}
+	// H diagonal over the eigen-transformed J, Java getWJtJlambda order (per
+	// parameter the contributions arrive in ascending value order; reg rows
+	// carry zero weight at this point and contribute exactly 0).
+	float h[POSE_NUM_PARAMS];
+	for (int par = 0; par < num_params; par++) h[par] = 0.0f;
+	for (int tile = 0; tile < num_tiles; tile++) {
+		const int vi = num_components * tile;
+		const float w = weights[vi];
+		if (!(w > 0.0f)) continue;
+		const bool have_raw = valid[tile] != 0;
+		const int ei = 4 * tile;
+		for (int par = 0; par < num_params; par++) {
+			const int raw_i = par * ns + vi;
+			const float jx = have_raw ? raw_jt[raw_i] : 0.0f;
+			const float jy = have_raw ? raw_jt[raw_i + 1] : 0.0f;
+			const float tj0 = pose_lma_add(pose_lma_mul(eigen[ei],     jx),
+			                               pose_lma_mul(eigen[ei + 1], jy));
+			const float tj1 = pose_lma_add(pose_lma_mul(eigen[ei + 2], jx),
+			                               pose_lma_mul(eigen[ei + 3], jy));
+			h[par] = pose_lma_add(h[par], pose_lma_mul(w, pose_lma_mul(tj0, tj0)));
+			h[par] = pose_lma_add(h[par], pose_lma_mul(w, pose_lma_mul(tj1, tj1)));
+		}
+	}
+	// Regularization weights (Java: param_regweights[i]/sqrt(WJtJ_ii), 0 if 0).
+	for (int par = 0; par < num_params; par++) {
+		const float rw = (h[par] > 0.0f) ?
+				__fdiv_rn(param_regweights[par], sqrtf(h[par])) : 0.0f;
+		weights[ns + par] = rw;
+		prep_result[POSE_LMA_PREP_REG_W0 + par] = rw;
+	}
+	// Full-weight normalization (Java normalizeWeights, serial ascending order).
+	float pure = 0.0f;
+	for (int i = 0; i < ns; i++) pure = pose_lma_add(pure, weights[i]);
+	float full = pure;
+	for (int par = 0; par < num_params; par++) full = pose_lma_add(full, weights[ns + par]);
+	const float s = __fdiv_rn(1.0f, full);
+	for (int i = 0; i < num_values; i++) weights[i] = pose_lma_mul(weights[i], s);
+	prep_result[POSE_LMA_PREP_PURE_WEIGHT] = __fdiv_rn(pure, full);
+	prep_result[POSE_LMA_PREP_FULL_WEIGHT] = full;
+	prep_result[POSE_LMA_PREP_VALID] = (isfinite(s) && (full > 0.0f)) ? 1.0f : 0.0f;
+}
+
 extern "C" __global__ void lma_normal_products(
 		int            num_params,
 		int            num_values,

src/tp_lma.h

@@ -39,6 +39,16 @@
 #define POSE_LMA_CANDIDATE_PURE_RMS 22
 #define POSE_LMA_ACCEPTED        23
 #define POSE_LMA_STOP            24
+// Resident prepare compact result (rung C1, design 2026-07-16). By Claude on 07/17/2026.
+#define POSE_LMA_PREP_RESULT      8
+#define POSE_LMA_PREP_VALID       0
+#define POSE_LMA_PREP_NUM_MEAS    1
+#define POSE_LMA_PREP_SUM_W_RAW   2
+#define POSE_LMA_PREP_PURE_WEIGHT 3
+#define POSE_LMA_PREP_FULL_WEIGHT 4
+#define POSE_LMA_PREP_REG_W0      5
+#define POSE_LMA_PREP_REG_W1      6
+#define POSE_LMA_PREP_REG_W2      7
 
 typedef struct {
 	int   width;
@@ -121,6 +131,67 @@ extern "C" __global__ void pose_lma_reduce_rms(
 		float       * result,         // candidate RMS fields in resident result
 		int         * current_set);   // device current/candidate set index, flipped on accept
 
+// ── Rung C1: resident prepareLMA (design 2026-07-16) ─────────────────────────
+// The measured inputs never leave the GPU: assemble scatters the resident peak
+// rows into the full-grid y/weights/eigen/selection buffers, then two SERIAL
+// one-thread tails replicate the Java float weight policy exactly (serial order
+// = bit-exact against the Java float-clone oracle; O(num_values) adds, ~us).
+// Kernel contracts stay parameterized (tile/corr counts, slot/shift, params);
+// CUAS specifics (150 tiles, 3 angles, slot 0xff) are host policy.
+// By Claude on 07/17/2026.
+
+// Scatter peak rows {dx,dy,str,e0x,e0y,l0,l1,valid} into per-tile LMA inputs.
+// y gets the MEASURED OFFSETS only (fx added by pose_lma_prepare_finish);
+// weights get the raw strength at [nc*tile] (component dup + normalization in
+// the tails); eigen gets the 2x2 row transform (setEigenTransform math, float)
+// or identity when use_eigen==0; selection marks measured tiles for the fx
+// kernel. Buffers must be pre-zeroed.
+extern "C" __global__ void pose_lma_assemble(
+		int                   num_corr,       // resident correlation rows
+		int                   num_tiles,      // full tile grid
+		int                   corr_slot,      // sum-slot id in the packed index (CUAS: 0xff)
+		int                   ntile_shift,    // packed-index tile shift (CUAS: 8)
+		int                   num_components, // 2 supported by this implementation
+		int                   use_eigen,      // 0 = identity transform for measured tiles
+		float                 eig_min_sqrt,
+		float                 eig_max_sqrt,
+		float                 min_confidence,
+		int                   same_weights,
+		const int           * corr_indices,   // [num_corr] packed (tile<<shift)|slot
+		const float         * peaks,          // [num_corr][8]
+		float               * y,              // [nc*num_tiles+num_params] pre-zeroed
+		float               * weights,        // [nc*num_tiles+num_params] pre-zeroed
+		float               * eigen,          // [4*num_tiles] pre-zeroed
+		unsigned char       * selection);     // [num_tiles] pre-zeroed
+
+// Serial tail 1: sample-weight normalization (sum 2w -> scale, duplicate the
+// second component) + pull-target y rows; counts measured tiles.
+extern "C" __global__ void pose_lma_prepare_norm(
+		int                   num_tiles,
+		int                   num_components,
+		int                   num_params,
+		const float         * pull_targets,   // [num_params]
+		float               * y,
+		float               * weights,
+		float               * prep_result);   // [POSE_LMA_PREP_RESULT]
+
+// Serial tail 2 (after pose_fx_jacobian at the current pose): y += fx for
+// measured+valid tiles; H diagonal from the eigen-transformed raw J (serial,
+// Java getWJtJlambda order); regularization weights; full-weight normalization
+// (pure_weight = sample fraction).
+extern "C" __global__ void pose_lma_prepare_finish(
+		int                   num_tiles,
+		int                   num_components,
+		int                   num_params,
+		const float         * raw_fx,
+		const float         * raw_jt,          // parameter-major [np][nc*num_tiles]
+		const unsigned char * valid,
+		const float         * eigen,
+		const float         * param_regweights,// [num_params]
+		float               * y,
+		float               * weights,
+		float               * prep_result);    // [POSE_LMA_PREP_RESULT]
+
 extern "C" __global__ void lma_normal_products(
 		int            num_params,
 		int            num_values,