CODEX: Continue pose LMA from resident Jacobians

Co-authored-by: Codex <codex@elphel.com>

src/main/java/com/elphel/imagej/cuas/rt/CuasPoseRT.java

@@ -836,7 +836,7 @@ public class CuasPoseRT {
 	private static boolean lma_products_path_reported = false;
 	// once-per-program console note for production getFxDerivs() on CUDA
 	private static boolean pose_fx_path_reported = false;
-	// once-per-program console note for the fixed float runLma candidate path
+	// once-per-program console note for resident float preparation/reduction/candidate
 	private static boolean pose_lma_step_path_reported = false;
 	// one raw-pixel/J comparison of the float oracle against the current double geometry
 	private static boolean float_jacobian_oracle_reported = false;
@@ -1034,16 +1034,19 @@ public class CuasPoseRT {
 			}
 			return result;
 		}); // first production gate: replace only getFxDerivs per-tile work; callers remain unchanged // By Codex on 07/15/2026
-		intersceneLma.setPoseLmaStepProvider((lambda, weights, jt, ymfxWeighted, vector) -> {
-			final IntersceneLmaFloat.LmaStepResult result = lmaGpu.execPoseLmaStep(
-					lambda, weights, jt, ymfxWeighted, vector);
+		intersceneLma.setPoseLmaResidentStepProvider((
+				lambda, weights, y, eigen, vector, pureWeight, numTiles, capturePrepared) -> {
+			final IntersceneLmaFloat.ResidentLmaStepResult result =
+					lmaGpu.execPoseLmaResidentStep(
+							lambda, weights, y, eigen, vector, pureWeight,
+							numTiles, capturePrepared);
 			if ((result != null) && !pose_lma_step_path_reported) {
 				pose_lma_step_path_reported = true;
-				System.out.println("CuasPoseRT: CUDA float runLma candidate active "+
-						"(one thread, per-call H2D/D2H; Java RMS/acceptance retained)");
+				System.out.println("CuasPoseRT: resident CUDA float LMA preparation/reduction active "+
+						"(resident raw fx/J reused; Java oracle readback and acceptance retained)");
 			}
 			return result;
-		}); // validation rung: only fixed 3x3 candidate calculation moves to CUDA // By Codex on 07/15/2026
+		}); // validation rung: resident raw fx/J -> parallel prepare/reduce -> fixed 3x3 tail // By Codex on 07/15/2026
 		intersceneLma.setNormalEquationProvider((weights, jt, ymfxWeighted) -> {
 			final long profileStart = (rtProfile != null) ? rtProfile.start() : 0L;
 			final double [] products = lmaGpu.execLmaNormalProducts(weights, jt, ymfxWeighted);

src/main/java/com/elphel/imagej/gpu/GpuQuad.java

@@ -2937,6 +2937,23 @@ public class GpuQuad{ // quad camera description
 		return null;
 	}
 
+	/**
+	 * Resident raw-fx/J continuation gate. JNA applies float eigen/residual
+	 * preparation and hierarchical products on-device, then runs the fixed 3x3
+	 * tail. Base/JCuda returns null so the Java solver remains the fallback.
+	 */
+	public IntersceneLmaFloat.ResidentLmaStepResult execPoseLmaResidentStep(
+			float lambda,
+			float [] weights,
+			float [] y,
+			float [] eigen,
+			float [] vector,
+			float pure_weight,
+			int num_tiles,
+			boolean capture_prepared) {
+		return null;
+	}
+
 	public void execRBGA(
 			double [] color_weights,
 			boolean   is_lwir,

src/main/java/com/elphel/imagej/gpu/jna/GpuQuadJna.java

@@ -553,6 +553,49 @@ public class GpuQuadJna extends GpuQuad {
         if (rc != 0) {
             throw new IllegalStateException("GpuQuadJna.execPoseLmaStep rc="+rc+": "+lib.tp_last_error());
         }
+        return unpackPoseLmaStep(packed);
+    }
+
+    // ---- resident float pose preparation/reduction continuation (roadmap rung 3-A4e) ----
+    @Override public IntersceneLmaFloat.ResidentLmaStepResult execPoseLmaResidentStep(
+            final float lambda,
+            final float[] weights,
+            final float[] y,
+            final float[] eigen,
+            final float[] vector,
+            final float pureWeight,
+            final int numTiles,
+            final boolean capturePrepared) {
+        final int numValues = IntersceneLmaFloat.NUM_COMPONENTS * numTiles +
+                IntersceneLmaFloat.NUM_PARAMS;
+        if ((numTiles <= 0) || (weights == null) || (weights.length != numValues) ||
+                (y == null) || (y.length != numValues) || (eigen == null) ||
+                (eigen.length != 4 * numTiles) || (vector == null) ||
+                (vector.length != IntersceneLmaFloat.NUM_PARAMS) ||
+                !Float.isFinite(pureWeight) || !(pureWeight > 0.0f)) {
+            throw new IllegalArgumentException("GpuQuadJna.execPoseLmaResidentStep: inconsistent inputs");
+        }
+        final float[] packed = new float[IntersceneLmaFloat.LMA_RESIDENT_RESULT_FLOATS];
+        final float[] preparedJt = capturePrepared ?
+                new float[IntersceneLmaFloat.NUM_PARAMS * numValues] : null;
+        final float[] preparedYmfx = capturePrepared ? new float[numValues] : null;
+        final int rc = lib.tp_proc_exec_pose_lma_resident_step(
+                proc, lambda, weights, y, eigen, vector, pureWeight, numTiles,
+                capturePrepared ? 1 : 0, packed, preparedJt, preparedYmfx);
+        if (rc != 0) {
+            throw new IllegalStateException("GpuQuadJna.execPoseLmaResidentStep rc="+rc+": "+
+                    lib.tp_last_error());
+        }
+        final IntersceneLmaFloat.PreparedLma prepared = capturePrepared ?
+                new IntersceneLmaFloat.PreparedLma(
+                        preparedJt, preparedYmfx,
+                        packed[IntersceneLmaFloat.LMA_RESULT_FLOATS],
+                        packed[IntersceneLmaFloat.LMA_RESULT_FLOATS + 1]) : null;
+        return new IntersceneLmaFloat.ResidentLmaStepResult(
+                unpackPoseLmaStep(packed), prepared);
+    }
+
+    private static IntersceneLmaFloat.LmaStepResult unpackPoseLmaStep(final float[] packed) {
         final int deltaOffset = 1 + IntersceneLmaFloat.LMA_PRODUCTS;
         final int candidateOffset = deltaOffset + IntersceneLmaFloat.NUM_PARAMS;
         return new IntersceneLmaFloat.LmaStepResult(

src/main/java/com/elphel/imagej/gpu/jna/Stage0.java

@@ -10,7 +10,7 @@ import com.sun.jna.Pointer;
  *        com.elphel.imagej.gpu.jna.Stage0 [kernel_src_dir] [libcudadevrt.a]
  */
 public class Stage0 {
-	private static final int EXPECTED_KERNELS = 30; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + pose step + products
+	private static final int EXPECTED_KERNELS = 32; // 19 base + 2 conditioning + 4 consolidate + 2 peak + pose fx/J + 3 pose LMA + products
     public static void main(String[] args) {
         String srcdir = (args.length > 0) ? args[0] : "/home/elphel/git/tile_processor_gpu/src";
         String devrt  = (args.length > 1) ? args[1] : "/usr/local/cuda/targets/x86_64-linux/lib/libcudadevrt.a";

src/main/java/com/elphel/imagej/gpu/jna/TpJna.java

@@ -12,8 +12,8 @@ public interface TpJna extends Library {
     /** NVRTC-compile the kernels in srcdir (+ getTpDefines), cuLink(libcudadevrt), load the module.
      *  Returns an opaque module handle, or null on failure (see tp_last_error()). */
     Pointer tp_create_module(String srcdir, String libcudadevrt);
-    /** Number of kernel functions resolved in the module (30 expected:
-     *  19 base + 2 conditioning + 4 consolidation + 2 peak + pose-fx/J + pose-step + LMA-products),
+    /** Number of kernel functions resolved in the module (32 expected:
+     *  19 base + 2 conditioning + 4 consolidation + 2 peak + pose-fx/J + 3 pose-LMA + LMA-products),
      *  or -1 if handle is null. */
     int     tp_module_num_functions(Pointer module);
     /** Last error message (NVRTC/CUDA log), empty if none. */
@@ -154,6 +154,14 @@ public interface TpJna extends Library {
                                       float lambda, float[] weights, float[] jt,
                                       float[] ymfxWeighted, float[] vector,
                                       int numValues, float[] result);
+    /** Continue from raw fx/J left by the most recent pose kernel. Result has
+     *  the 19 candidate floats followed by current full/pure RMS. Prepared
+     *  arrays are nullable unless capturePrepared is nonzero. */
+    int     tp_proc_exec_pose_lma_resident_step(Pointer proc,
+                                      float lambda, float[] weights, float[] y,
+                                      float[] eigen, float[] vector, float pureWeight,
+                                      int numTiles, int capturePrepared, float[] result,
+                                      float[] preparedJt, float[] preparedYmfx);
     /** Deterministic double normal-equation products for the lean pose LMA.
      *  jt is parameter-major [numParams][numValues]; out is row-major H followed by b.
      *  Damping and the small solve remain in Java. 0 on success. */

src/main/java/com/elphel/imagej/tileprocessor/IntersceneLma.java

@@ -111,11 +111,26 @@ public class IntersceneLma {
 				float [] vector);
 	}
 	private PoseLmaStepProvider poseLmaStepProvider = null;
+	/** Optional continuation from raw fx/J already resident in the pose backend. */
+	@FunctionalInterface
+	public interface PoseLmaResidentStepProvider {
+		IntersceneLmaFloat.ResidentLmaStepResult calculate(
+				float lambda,
+				float [] weights,
+				float [] y,
+				float [] eigen,
+				float [] vector,
+				float pureWeight,
+				int numTiles,
+				boolean capturePrepared);
+	}
+	private PoseLmaResidentStepProvider poseLmaResidentStepProvider = null;
 	private static boolean poseLmaStepOracleReported = false;
 	private IntersceneLmaFloat.Camera poseFloatReference = null;
 	private IntersceneLmaFloat.Camera poseFloatScene = null;
 	private float [] poseFloatCenters = null;
 	private boolean [] poseFloatSelection = null;
+	private IntersceneLmaFloat.Result poseFloatLastResult = null;
 	
 	public IntersceneLma(
 			boolean thread_invariant,
@@ -141,11 +156,16 @@ public class IntersceneLma {
 	public void setPoseLmaStepProvider(final PoseLmaStepProvider provider) {
 		this.poseLmaStepProvider = provider;
 	}
+	/** Install the resident raw-fx/J continuation backend. */
+	public void setPoseLmaResidentStepProvider(final PoseLmaResidentStepProvider provider) {
+		this.poseLmaResidentStepProvider = provider;
+	}
 	private void resetPoseFxInputs() {
 		poseFloatReference = null;
 		poseFloatScene = null;
 		poseFloatCenters = null;
 		poseFloatSelection = null;
+		poseFloatLastResult = null;
 	}
 	public double [][]       getLastJT(){
 		return last_jt;
@@ -976,38 +996,62 @@ public class IntersceneLma {
 		for (int i = 0; i < parameters_vector.length; i++) {
 			new_vector[i] += scale * delta[i];
 		}
-		if ((poseLmaStepProvider != null) && isPoseLmaStepShape()) {
+		if (((poseLmaResidentStepProvider != null) || (poseLmaStepProvider != null)) &&
+				isPoseLmaStepShape()) {
 			final float [] floatWeights = toFloatArray(weights);
-			final float [] floatJt = flattenJtFloat(last_jt);
-			final float [] floatYmfx = toFloatArray(last_ymfx);
 			final float [] floatVector = toFloatArray(parameters_vector);
-			final IntersceneLmaFloat.LmaStepResult poseStep = poseLmaStepProvider.calculate(
-					(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
-			if ((debug_level >= 1) && !poseLmaStepOracleReported && (poseStep != null)) {
-				poseLmaStepOracleReported = true;
-				final IntersceneLmaFloat.LmaStepResult javaFloatStep = IntersceneLmaFloat.lmaStep(
-						(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
-				final IntersceneLmaFloat.LmaStepComparison comparison =
-						IntersceneLmaFloat.compareLmaSteps(poseStep, javaFloatStep);
-				System.out.println("IntersceneLma CUDA vs Java-float runLma candidate: "+
-						comparison.format());
-				float maxProducts = 0.0f;
-				if ((normalProducts != null) && (normalProducts.length == poseStep.products.length)) {
-					for (int i = 0; i < normalProducts.length; i++) {
-						maxProducts = Math.max(maxProducts,
-								Math.abs(poseStep.products[i] - (float) normalProducts[i]));
+			IntersceneLmaFloat.LmaStepResult poseStep = null;
+			if ((poseLmaResidentStepProvider != null) && (poseFloatLastResult != null) &&
+					(poseFloatLastResult.jt != null)) {
+				final boolean capturePrepared = (debug_level >= 1) && !poseLmaStepOracleReported;
+				final float [] floatY = toFloatArray(y_vector);
+				final float [] floatEigen = flattenEigenTransform(
+						eig_trans, poseFloatLastResult.numTiles);
+				final IntersceneLmaFloat.ResidentLmaStepResult resident =
+						poseLmaResidentStepProvider.calculate(
+								(float) lambda, floatWeights, floatY, floatEigen, floatVector,
+								(float) pure_weight, poseFloatLastResult.numTiles, capturePrepared);
+				if (resident != null) poseStep = resident.step;
+				if (capturePrepared && (resident != null)) {
+					if (resident.prepared == null) {
+						throw new IllegalStateException("IntersceneLma: resident LMA oracle arrays missing");
 					}
+					final IntersceneLmaFloat.PreparedLma javaPrepared =
+							IntersceneLmaFloat.prepareLma(
+									poseFloatLastResult, floatY, floatWeights, floatEigen,
+									floatVector, (float) pure_weight);
+					final IntersceneLmaFloat.PreparationComparison preparationComparison =
+							IntersceneLmaFloat.comparePreparation(resident.prepared, javaPrepared);
+					System.out.println("IntersceneLma resident CUDA vs Java-float preparation: "+
+							preparationComparison.format());
+					final IntersceneLmaFloat.LmaStepResult javaFloatStep = IntersceneLmaFloat.lmaStep(
+							(float) lambda, floatWeights, javaPrepared.jt,
+							javaPrepared.ymfxWeighted, floatVector);
+					final IntersceneLmaFloat.LmaStepComparison comparison =
+							IntersceneLmaFloat.compareLmaSteps(poseStep, javaFloatStep);
+					System.out.println("IntersceneLma resident hierarchical vs Java-float serial candidate: "+
+							comparison.format());
+					reportPoseStepVsDouble("resident CUDA-float", poseStep,
+							normalProducts, delta, new_vector);
+					poseLmaStepOracleReported = true;
 				}
-				float maxDelta = 0.0f;
-				float maxCandidate = 0.0f;
-				for (int i = 0; i < IntersceneLmaFloat.NUM_PARAMS; i++) {
-					maxDelta = Math.max(maxDelta, Math.abs(poseStep.delta[i] - (float) delta[i]));
-					maxCandidate = Math.max(maxCandidate,
-							Math.abs(poseStep.candidate[i] - (float) new_vector[i]));
+			}
+			if ((poseStep == null) && (poseLmaStepProvider != null)) {
+				final float [] floatJt = flattenJtFloat(last_jt);
+				final float [] floatYmfx = toFloatArray(last_ymfx);
+				poseStep = poseLmaStepProvider.calculate(
+						(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
+				if ((debug_level >= 1) && !poseLmaStepOracleReported && (poseStep != null)) {
+					final IntersceneLmaFloat.LmaStepResult javaFloatStep = IntersceneLmaFloat.lmaStep(
+							(float) lambda, floatWeights, floatJt, floatYmfx, floatVector);
+					final IntersceneLmaFloat.LmaStepComparison comparison =
+							IntersceneLmaFloat.compareLmaSteps(poseStep, javaFloatStep);
+					System.out.println("IntersceneLma CUDA vs Java-float runLma candidate: "+
+							comparison.format());
+					reportPoseStepVsDouble("CUDA-float", poseStep,
+							normalProducts, delta, new_vector);
+					poseLmaStepOracleReported = true;
 				}
-				System.out.println(String.format(
-						"IntersceneLma CUDA-float vs Java-double runLma candidate: valid=%s, max|H,b|=%g, max|delta|=%g, max|candidate|=%g",
-						poseStep.valid, maxProducts, maxDelta, maxCandidate));
 			}
 			if ((poseStep != null) && poseStep.valid &&
 					(poseStep.candidate.length == new_vector.length)) {
@@ -1811,6 +1855,7 @@ public class IntersceneLma {
 						(result.jt.length != IntersceneLmaFloat.NUM_PARAMS * 2 * result.numTiles)))) {
 			throw new IllegalStateException("IntersceneLma: pose fx/J backend returned inconsistent arrays");
 		}
+		poseFloatLastResult = result;
 		for (int tile = 0; tile < result.numTiles; tile++) {
 			if (result.valid[tile] == 0) continue; // legacy geometry leaves zero rows for an unusable LMA tile
 			final int vi = 2 * tile;
@@ -1863,6 +1908,60 @@ public class IntersceneLma {
 		return result;
 	}
 
+	private static float [] flattenEigenTransform(
+			final double [][][] source,
+			final int numTiles) {
+		final float [] result = new float [4 * numTiles];
+		if (source == null) {
+			for (int tile = 0; tile < numTiles; tile++) {
+				result[4 * tile] = 1.0f;
+				result[4 * tile + 3] = 1.0f;
+			}
+			return result;
+		}
+		if (source.length != numTiles) {
+			throw new IllegalArgumentException("IntersceneLma: eigen transform tile count");
+		}
+		for (int tile = 0; tile < numTiles; tile++) {
+			if (source[tile] == null) continue;
+			if ((source[tile].length != 2) || (source[tile][0] == null) ||
+					(source[tile][1] == null) || (source[tile][0].length != 2) ||
+					(source[tile][1].length != 2)) {
+				throw new IllegalArgumentException("IntersceneLma: eigen transform["+tile+"]");
+			}
+			result[4 * tile] = (float) source[tile][0][0];
+			result[4 * tile + 1] = (float) source[tile][0][1];
+			result[4 * tile + 2] = (float) source[tile][1][0];
+			result[4 * tile + 3] = (float) source[tile][1][1];
+		}
+		return result;
+	}
+
+	private static void reportPoseStepVsDouble(
+			final String label,
+			final IntersceneLmaFloat.LmaStepResult poseStep,
+			final double [] normalProducts,
+			final double [] delta,
+			final double [] candidate) {
+		float maxProducts = 0.0f;
+		if ((normalProducts != null) && (normalProducts.length == poseStep.products.length)) {
+			for (int i = 0; i < normalProducts.length; i++) {
+				maxProducts = Math.max(maxProducts,
+						Math.abs(poseStep.products[i] - (float) normalProducts[i]));
+			}
+		}
+		float maxDelta = 0.0f;
+		float maxCandidate = 0.0f;
+		for (int i = 0; i < IntersceneLmaFloat.NUM_PARAMS; i++) {
+			maxDelta = Math.max(maxDelta, Math.abs(poseStep.delta[i] - (float) delta[i]));
+			maxCandidate = Math.max(maxCandidate,
+					Math.abs(poseStep.candidate[i] - (float) candidate[i]));
+		}
+		System.out.println(String.format(
+				"IntersceneLma %s vs Java-double runLma candidate: valid=%s, max|H,b|=%g, max|delta|=%g, max|candidate|=%g",
+				label, poseStep.valid, maxProducts, maxDelta, maxCandidate));
+	}
+
 	private double [][] getWJtJlambda( // USED in lwir
 			final double      lambda,
 			final double [][] jt)

src/main/java/com/elphel/imagej/tileprocessor/IntersceneLmaFloat.java

@@ -42,6 +42,7 @@ public final class IntersceneLmaFloat {
 	public static final int NUM_COMPONENTS = 2;
 	public static final int LMA_PRODUCTS = NUM_PARAMS * NUM_PARAMS + NUM_PARAMS;
 	public static final int LMA_RESULT_FLOATS = 1 + LMA_PRODUCTS + 2 * NUM_PARAMS;
+	public static final int LMA_RESIDENT_RESULT_FLOATS = LMA_RESULT_FLOATS + 2;
 	public static final float INFINITY_DISPARITY = 0.01f;
 
 	static final int ERS_XYZ = 0;
@@ -192,6 +193,51 @@ public final class IntersceneLmaFloat {
 		}
 	}
 
+	/** Float eigen/residual preparation immediately downstream of raw resident fx/J. */
+	public static final class PreparedLma {
+		public final float [] jt;
+		public final float [] ymfxWeighted;
+		public final float rms;
+		public final float pureRms;
+
+		public PreparedLma(
+				final float [] jt,
+				final float [] ymfxWeighted,
+				final float rms,
+				final float pureRms) {
+			this.jt = jt;
+			this.ymfxWeighted = ymfxWeighted;
+			this.rms = rms;
+			this.pureRms = pureRms;
+		}
+	}
+
+	/** Compact resident step plus optional prepared arrays for a one-shot oracle gate. */
+	public static final class ResidentLmaStepResult {
+		public final LmaStepResult step;
+		public final PreparedLma prepared;
+
+		public ResidentLmaStepResult(final LmaStepResult step, final PreparedLma prepared) {
+			this.step = step;
+			this.prepared = prepared;
+		}
+	}
+
+	/** Statistics for the resident preparation arrays and current-point RMS. */
+	public static final class PreparationComparison {
+		public boolean bitExact;
+		public float maxAbsJt;
+		public float maxAbsYmfx;
+		public float absRms;
+		public float absPureRms;
+
+		public String format() {
+			return String.format(
+					"bit-exact=%s, max|J|=%g, max|W(y-f)|=%g, |RMS|=%g, |pure RMS|=%g",
+					bitExact, maxAbsJt, maxAbsYmfx, absRms, absPureRms);
+		}
+	}
+
 	/**
 	 * Primitive-float clone of the fixed three-angle candidate calculation in
 	 * {@code IntersceneLma.lmaStep()}. It deliberately keeps all reductions and
@@ -249,6 +295,75 @@ public final class IntersceneLmaFloat {
 		return new LmaStepResult(allFinite(delta) && allFinite(candidate), products, delta, candidate);
 	}
 
+	/**
+	 * Primitive-float oracle for the resident per-tile preparation kernel. Raw
+	 * projection/Jacobian values are eigen-transformed, weighted residuals are
+	 * formed, and the three pull rows are appended without any general matrix API.
+	 */
+	public static PreparedLma prepareLma(
+			final Result raw,
+			final float [] y,
+			final float [] weights,
+			final float [] eigen,
+			final float [] vector,
+			final float pureWeight) {
+		if ((raw == null) || (raw.numTiles <= 0) || (raw.fx == null) ||
+				(raw.fx.length != NUM_COMPONENTS * raw.numTiles) || (raw.jt == null) ||
+				(raw.jt.length != NUM_PARAMS * NUM_COMPONENTS * raw.numTiles) ||
+				(raw.valid == null) || (raw.valid.length != raw.numTiles) ||
+				(y == null) || (weights == null) ||
+				(y.length != NUM_COMPONENTS * raw.numTiles + NUM_PARAMS) ||
+				(weights.length != y.length) || (eigen == null) ||
+				(eigen.length != 4 * raw.numTiles) || (vector == null) ||
+				(vector.length != NUM_PARAMS) || !finite(pureWeight) || !(pureWeight > 0.0f)) {
+			throw new IllegalArgumentException("inconsistent resident pose LMA preparation arrays");
+		}
+		final int numValues = y.length;
+		final float [] transformedJt = new float [NUM_PARAMS * numValues];
+		final float [] ymfx = new float [numValues];
+		float pureL2 = 0.0f;
+		for (int tile = 0; tile < raw.numTiles; tile++) {
+			final int vi = NUM_COMPONENTS * tile;
+			final int ei = 4 * tile;
+			final boolean useTile = weights[vi] > 0.0f;
+			final boolean haveRaw = useTile && (raw.valid[tile] != 0);
+			final float rawX = haveRaw ? raw.fx[vi] : 0.0f;
+			final float rawY = haveRaw ? raw.fx[vi + 1] : 0.0f;
+			float residual0 = 0.0f;
+			float residual1 = 0.0f;
+			if (useTile) {
+				final float dx = y[vi] - rawX;
+				final float dy = y[vi + 1] - rawY;
+				residual0 = eigen[ei] * dx + eigen[ei + 1] * dy;
+				residual1 = eigen[ei + 2] * dx + eigen[ei + 3] * dy;
+			}
+			ymfx[vi] = residual0 * weights[vi];
+			ymfx[vi + 1] = residual1 * weights[vi + 1];
+			pureL2 += residual0 * ymfx[vi];
+			pureL2 += residual1 * ymfx[vi + 1];
+			for (int par = 0; par < NUM_PARAMS; par++) {
+				final int rawIndex = par * NUM_COMPONENTS * raw.numTiles + vi;
+				final float jx = haveRaw ? raw.jt[rawIndex] : 0.0f;
+				final float jy = haveRaw ? raw.jt[rawIndex + 1] : 0.0f;
+				final int outIndex = par * numValues + vi;
+				transformedJt[outIndex] = eigen[ei] * jx + eigen[ei + 1] * jy;
+				transformedJt[outIndex + 1] = eigen[ei + 2] * jx + eigen[ei + 3] * jy;
+			}
+		}
+		float fullL2 = pureL2;
+		for (int par = 0; par < NUM_PARAMS; par++) {
+			final int value = NUM_COMPONENTS * raw.numTiles + par;
+			transformedJt[par * numValues + value] = 1.0f;
+			final float residual = y[value] - vector[par];
+			ymfx[value] = residual * weights[value];
+			fullL2 += residual * ymfx[value];
+		}
+		return new PreparedLma(
+				transformedJt, ymfx,
+				(float) Math.sqrt(fullL2),
+				(float) Math.sqrt(pureL2 / pureWeight));
+	}
+
 	public static LmaStepComparison compareLmaSteps(
 			final LmaStepResult actual,
 			final LmaStepResult expected) {
@@ -270,6 +385,43 @@ public final class IntersceneLmaFloat {
 		return comparison;
 	}
 
+	public static PreparationComparison comparePreparation(
+			final PreparedLma actual,
+			final PreparedLma expected) {
+		if ((actual == null) || (expected == null) || (actual.jt == null) ||
+				(expected.jt == null) || (actual.jt.length != expected.jt.length) ||
+				(actual.ymfxWeighted == null) || (expected.ymfxWeighted == null) ||
+				(actual.ymfxWeighted.length != expected.ymfxWeighted.length)) {
+			throw new IllegalArgumentException("incompatible resident pose LMA preparation results");
+		}
+		final PreparationComparison comparison = new PreparationComparison();
+		comparison.bitExact = true;
+		comparison.maxAbsJt = comparePreparedArrays(actual.jt, expected.jt, comparison);
+		comparison.maxAbsYmfx = comparePreparedArrays(
+				actual.ymfxWeighted, expected.ymfxWeighted, comparison);
+		comparison.absRms = Math.abs(actual.rms - expected.rms);
+		comparison.absPureRms = Math.abs(actual.pureRms - expected.pureRms);
+		return comparison;
+	}
+
+	private static float comparePreparedArrays(
+			final float [] actual,
+			final float [] expected,
+			final PreparationComparison comparison) {
+		float maximum = 0.0f;
+		for (int i = 0; i < actual.length; i++) {
+			final boolean sameBits = Float.floatToRawIntBits(actual[i]) ==
+					Float.floatToRawIntBits(expected[i]);
+			if (!sameBits) comparison.bitExact = false;
+			if (finite(actual[i]) && finite(expected[i])) {
+				maximum = Math.max(maximum, Math.abs(actual[i] - expected[i]));
+			} else if (!sameBits) {
+				maximum = Float.POSITIVE_INFINITY;
+			}
+		}
+		return maximum;
+	}
+
 	private static float compareFloatArrays(
 			final float [] actual,
 			final float [] expected,