CODEX: Move lean pose normal products to GPU

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

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

@@ -525,6 +525,8 @@ public class CuasPoseRT {
 
 	// once-per-program console note that the GPU peak path is live // By Claude on 07/14/2026
 	private static boolean peaks_path_reported = false;
+	// once-per-program console note for the rung-3 normal-equation product path
+	private static boolean lma_products_path_reported = false;
 
 	private static int GPUTileProcessorDttSize() {
 		return com.elphel.imagej.gpu.GPUTileProcessor.DTT_SIZE;
@@ -702,6 +704,16 @@ public class CuasPoseRT {
 		final IntersceneLma intersceneLma = new IntersceneLma(
 				clt_parameters.ilp.ilma_thread_invariant,
 				0.0); // no disparity weight (2D only)
+		final GpuQuad lmaGpu = center_CLT.getGPUQuad();
+		intersceneLma.setNormalEquationProvider((weights, jt, ymfxWeighted) -> {
+			final double [] products = lmaGpu.execLmaNormalProducts(weights, jt, ymfxWeighted);
+			if ((products != null) && !lma_products_path_reported) {
+				lma_products_path_reported = true;
+				System.out.println("CuasPoseRT: GPU LMA normal-equation products active "+
+						"(lma_normal_products, double) - 3x3 solve stays on CPU");
+			}
+			return products;
+		}); // roadmap rung 3: J^T W J + J^T W(y-f) on GPU; solve/acceptance stay CPU // By Codex on 07/14/2026
 		double [][] scene_xyzatr0 = new double [][] {predicted[0].clone(), predicted[1].clone()};
 		double [][][] cm = null;
 		// lean v2 (Andrey's ruling 07/12/2026): pose_cycles > 0 = run EXACTLY that many

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

@@ -2812,6 +2812,21 @@ public class GpuQuad{ // quad camera description
 		return false;
 	}
 
+	// ---- LMA normal-equation products (tp_lma.cu, roadmap rung 3) ----
+	/**
+	 * Compute H=J^T W J and b=J^T W(y-f) on backends that carry the rung-3
+	 * kernel. jt is parameter-major; ymfx_weighted already includes W and any
+	 * eigen transform. The returned flat array is row-major H followed by b.
+	 * Damping, inversion, candidate evaluation and acceptance remain on CPU.
+	 * Base/JCuda returns null (unsupported).
+	 */
+	public double [] execLmaNormalProducts(
+			double [] weights,
+			double [][] jt,
+			double [] ymfx_weighted) {
+		return null;
+	}
+
 	public void execRBGA(
 			double [] color_weights,
 			boolean   is_lwir,
@@ -5360,4 +5375,4 @@ public class GpuQuad{ // quad camera description
 	}
 
 
-} // end of public class GpuQuad
\ No newline at end of file
+} // end of public class GpuQuad

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

@@ -80,6 +80,11 @@ public class GpuQuadJna extends GpuQuad {
     // clobbering any explicitly uploaded data (setBayerImages(data, true) / setBayerImage()) - found by the
     // curt_cond_test raw baseline coming out bit-identical to the conditioned render.
     private boolean jna_bayer_set = false;
+    // Native equivalent of GpuQuad.geometry_correction_set. The reverse-distortion
+    // table depends only on GeometryCorrection, not on the per-iteration pose/tasks;
+    // keep the upload/table cached until resetGeometryCorrection() invalidates it.
+    // By Codex on 07/14/2026 (roadmap rung-3 hoist).
+    private boolean jna_geometry_set = false;
 
     // Switch this shared GpuQuad to a different scene's QuadCLT. Base also clears
     // gpuTileProcessor.bayer_set (null here) — mirrored with jna_bayer_set: the next convert re-uploads.
@@ -90,6 +95,10 @@ public class GpuQuadJna extends GpuQuad {
         resetGeometryCorrectionVector();
         jna_bayer_set = false; // By Claude on 07/02/2026: match base (scene switch invalidates uploaded bayer)
     }
+    @Override public void resetGeometryCorrection() {
+        super.resetGeometryCorrection();
+        jna_geometry_set = false;
+    }
     @Override public void resetBayer() { jna_bayer_set = false; } // By Claude on 07/02/2026: match base semantics
 
     // CLT sizing (base derefs null gpu_clt_wh) — full-frame dims; matches TpProc's slice.
@@ -190,10 +199,16 @@ public class GpuQuadJna extends GpuQuad {
     }
 
     // ---- geometry ----
-    @Override public void setGeometryCorrection() { setGeometryCorrection(quadCLT.getGeometryCorrection(), false); }
+    @Override public void setGeometryCorrection() {
+        if (!jna_geometry_set) setGeometryCorrection(quadCLT.getGeometryCorrection(), false);
+    }
     @Override public void setGeometryCorrection(GeometryCorrection gc, boolean use_java_rByRDist) {
         float[] fgc = gc.expandSensors(GPUTileProcessor.MAX_NUM_CAMS).toFloatArray();
-        lib.tp_proc_set_geometry(proc, fgc, fgc.length);
+        int rc = lib.tp_proc_set_geometry(proc, fgc, fgc.length);
+        if (rc != 0) {
+            throw new IllegalStateException("GpuQuadJna.setGeometryCorrection rc="+rc+": "+lib.tp_last_error());
+        }
+        jna_geometry_set = true;
     }
     @Override public void setExtrinsicsVector(CorrVector cv) {
         double[] dcv = cv.toFullRollArray();
@@ -382,6 +397,31 @@ public class GpuQuadJna extends GpuQuad {
         return lib.tp_proc_set_peak_debias(proc, debias, (debias == null) ? 0 : debias.length) == 0;
     }
 
+    // ---- LMA normal-equation products (tp_lma.cu, roadmap rung 3) ---- // By Codex on 07/14/2026
+    @Override public double[] execLmaNormalProducts(
+            double[] weights, double[][] jt, double[] ymfxWeighted) {
+        if ((weights == null) || (jt == null) || (jt.length == 0) ||
+                (ymfxWeighted == null) || (ymfxWeighted.length != weights.length)) {
+            throw new IllegalArgumentException("GpuQuadJna.execLmaNormalProducts: inconsistent inputs");
+        }
+        final int numParams = jt.length;
+        final int numValues = weights.length;
+        final double[] flatJt = new double[numParams * numValues];
+        for (int i = 0; i < numParams; i++) {
+            if ((jt[i] == null) || (jt[i].length != numValues)) {
+                throw new IllegalArgumentException("GpuQuadJna.execLmaNormalProducts: jt["+i+"] length");
+            }
+            System.arraycopy(jt[i], 0, flatJt, i * numValues, numValues);
+        }
+        final double[] out = new double[numParams * numParams + numParams];
+        final int rc = lib.tp_proc_exec_lma_products(
+                proc, weights, flatJt, ymfxWeighted, numParams, numValues, out);
+        if (rc != 0) {
+            throw new IllegalStateException("GpuQuadJna.execLmaNormalProducts rc="+rc+": "+lib.tp_last_error());
+        }
+        return out;
+    }
+
     @Override public float[][] getCorr2D(int corr_rad) {
         int corr_size = (2 * corr_rad + 1) * (2 * corr_rad + 1);
         int n = lib.tp_proc_num_corr_tiles(proc);

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

+/**
+ ** LmaProducts.java - JNA/Java gate for the rung-3 LMA normal-equation kernel
+ **
+ ** Copyright (C) 2026 Elphel, Inc.
+ **
+ ** -----------------------------------------------------------------------------**
+ **
+ **  LmaProducts.java is free software: you can redistribute it and/or modify
+ **  it under the terms of the GNU General Public License as published by
+ **  the Free Software Foundation, either version 3 of the License, or
+ **  (at your option) any later version.
+ **
+ **  LmaProducts.java is distributed in the hope that it will be useful,
+ **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ **  GNU General Public License for more details.
+ **
+ **  You should have received a copy of the GNU General Public License
+ **  along with LmaProducts.java.  If not, see <http://www.gnu.org/licenses/>.
+ ** -----------------------------------------------------------------------------**
+ **
+ */
+package com.elphel.imagej.gpu.jna;
+
+/**
+ * Full Java -> JNA -> NVRTC gate for {@code lma_normal_products}. It uses the
+ * current lean shape (3 parameters, 2*150 samples + 3 regularization rows) and
+ * compares every returned double bit-for-bit with Java's own expression order.
+ *
+ * Run after building libtileproc.so:
+ * <pre>
+ * java -Djna.library.path=/home/elphel/git/tile_processor_gpu/jna \
+ *   -cp target/classes:~/.m2/repository/net/java/dev/jna/jna/5.14.0/jna-5.14.0.jar \
+ *   com.elphel.imagej.gpu.jna.LmaProducts \
+ *   /home/elphel/git/tile_processor_gpu/src \
+ *   /usr/local/cuda-12.8/targets/x86_64-linux/lib/libcudadevrt.a
+ * </pre>
+ * By Codex on 07/14/2026.
+ */
+public class LmaProducts {
+	public static void main(final String [] args) {
+		final String srcdir = (args.length > 0) ? args[0] : "/home/elphel/git/tile_processor_gpu/src";
+		final String devrt = (args.length > 1) ? args[1] :
+				"/usr/local/cuda-12.8/targets/x86_64-linux/lib/libcudadevrt.a";
+		final int numParams = 3;
+		final int numValues = 303;
+		final double [] weights = new double [numValues];
+		final double [][] jt = new double [numParams][numValues];
+		final double [] ymfx = new double [numValues];
+		for (int k = 0; k < numValues; k++) {
+			weights[k] = ((k % 11) == 0) ? 0.0 : (0.00037 * (k + 1) / (k + 3.25));
+			ymfx[k] = ((k & 1) != 0 ? -1.0 : 1.0) * (0.017 + 0.000031 * k * k);
+			for (int i = 0; i < numParams; i++) {
+				final double sign = (((k + 2 * i) & 1) != 0) ? -1.0 : 1.0;
+				jt[i][k] = sign * (13.7 * (i + 1) + 0.00091 * (k + 1) +
+						0.0000013 * ((k * (i + 3)) % 17));
+			}
+		}
+		final double [] expected = new double [numParams * numParams + numParams];
+		for (int i = 0; i < numParams; i++) {
+			for (int j = i; j < numParams; j++) {
+				double sum = 0.0;
+				for (int k = 0; k < numValues; k++) {
+					sum += weights[k] * jt[i][k] * jt[j][k];
+				}
+				expected[i * numParams + j] = expected[j * numParams + i] = sum;
+			}
+			double sum = 0.0;
+			for (int k = 0; k < numValues; k++) {
+				sum += jt[i][k] * ymfx[k];
+			}
+			expected[numParams * numParams + i] = sum;
+		}
+
+		int bad = 0;
+		final GpuQuadJna gpu = new GpuQuadJna(16, 16, 1, srcdir, devrt, -1);
+		try {
+			final double [] got = gpu.execLmaNormalProducts(weights, jt, ymfx);
+			for (int i = 0; i < got.length; i++) {
+				if (Double.doubleToRawLongBits(got[i]) != Double.doubleToRawLongBits(expected[i])) {
+					System.out.println("["+i+"] got "+got[i]+" expected "+expected[i]+" -> FAIL");
+					bad++;
+				}
+			}
+		} finally {
+			gpu.close();
+		}
+		System.out.println("Java/JNA lma_normal_products: "+expected.length+
+				" doubles, bit-exact -> "+((bad == 0) ? "PASS" : "FAIL"));
+		if (bad != 0) System.exit(1);
+	}
+}

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

@@ -4,12 +4,13 @@ import com.sun.jna.Native;
 import com.sun.jna.Pointer;
 
 /**
- * Stage-0b proof: from Java, load libtileproc.so via JNA, call tp_create_module, assert 19 kernels.
+ * Stage-0b proof: from Java, load libtileproc.so via JNA, call tp_create_module, assert all kernels.
  * Same result as the C++ probe but across the JNA boundary on the RTX 5060 Ti. By Claude on 2026-06-25.
  *   java -Djna.library.path=<dir with libtileproc.so> -cp target/classes:<jna.jar> \
  *        com.elphel.imagej.gpu.jna.Stage0 [kernel_src_dir] [libcudadevrt.a]
  */
 public class Stage0 {
+	private static final int EXPECTED_KERNELS = 26; // 19 base + 4 consolidate + 2 peak + 1 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";
@@ -24,7 +25,9 @@ public class Stage0 {
         String err = lib.tp_last_error();
         System.out.println("Stage0 (JNA): module created, functions loaded = " + n + (err != null && !err.isEmpty() ? ("  [warn: " + err + "]") : ""));
         lib.tp_destroy_module(m);
-        System.out.println(n == 19 ? "RESULT: PASS (19/19 via JNA)" : ("RESULT: FAIL (" + n + "/19)"));
-        System.exit(n == 19 ? 0 : 2);
+        System.out.println(n == EXPECTED_KERNELS ?
+				"RESULT: PASS ("+EXPECTED_KERNELS+"/"+EXPECTED_KERNELS+" via JNA)" :
+				("RESULT: FAIL (" + n + "/"+EXPECTED_KERNELS+")"));
+        System.exit(n == EXPECTED_KERNELS ? 0 : 2);
     }
 }

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 (23 expected: 19 + the 4 tp_consolidate
-     *  kernels added 07/12/2026), or -1 if handle is null. */
+    /** Number of kernel functions resolved in the module (26 expected through rung 3:
+     *  19 base + 4 consolidation + 2 peak + 1 LMA-products), or -1 if handle is null. */
     int     tp_module_num_functions(Pointer module);
     /** Last error message (NVRTC/CUDA log), empty if none. */
     String  tp_last_error();
@@ -120,6 +120,12 @@ public interface TpJna extends Library {
                                      int useFloat);
     /** D2H the peak rows (count x 8 floats). Returns the count. */
     int     tp_proc_get_peaks(Pointer proc, float[] out);
+    /** 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. */
+    int     tp_proc_exec_lma_products(Pointer proc,
+                                      double[] weights, double[] jt, double[] ymfxWeighted,
+                                      int numParams, int numValues, double[] out);
     int     tp_proc_num_corr_tiles(Pointer proc);
     int     tp_proc_num_corr_combo(Pointer proc);
     /** Upload to a named __constant__ symbol (lpf_data / lpf_corr / lpf_rb_corr). 0 on success. */

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

@@ -78,6 +78,13 @@ public class IntersceneLma {
 	private double [][][]     eig_trans = null;
 	private int               tilesX = -1;
 	private String            dbg_prefix = null;
+	/** Optional backend for the two normal-equation products in lmaStep().
+	 *  Inputs retain the Java/oracle representation; output is row-major H then b. */
+	@FunctionalInterface
+	public interface NormalEquationProvider {
+		double [] products(double [] weights, double [][] jt, double [] ymfxWeighted);
+	}
+	private NormalEquationProvider normalEquationProvider = null;
 	
 	public IntersceneLma(
 			boolean thread_invariant,
@@ -91,6 +98,10 @@ public class IntersceneLma {
 	public int getNumComponents() {
 		return num_components;
 	}
+	/** Install a normal-equation product backend. Null restores the CPU implementation. */
+	public void setNormalEquationProvider(final NormalEquationProvider provider) {
+		this.normalEquationProvider = provider;
+	}
 	public double [][]       getLastJT(){
 		return last_jt;
 	}
@@ -850,11 +861,33 @@ public class IntersceneLma {
 				showDebugImage(dbg_prefix+"-INIT");
 			}
 		}
-		Matrix y_minus_fx_weighted = new Matrix(this.last_ymfx, this.last_ymfx.length);
-
-		Matrix wjtjlambda = new Matrix(getWJtJlambda(
-				lambda, // *10, // temporary
-				this.last_jt)); // double [][] jt)
+		Matrix wjtjlambda;
+		Matrix jty;
+		if (normalEquationProvider != null) {
+			final int numPars = this.last_jt.length;
+			final double [] products = normalEquationProvider.products(
+					this.weights, this.last_jt, this.last_ymfx);
+			if ((products == null) || (products.length != numPars * numPars + numPars)) {
+				throw new IllegalStateException("IntersceneLma: normal-equation provider returned "+
+						((products == null) ? "null" : products.length)+" values, expected "+
+						(numPars * numPars + numPars));
+			}
+			final double [][] h = new double [numPars][numPars];
+			for (int i = 0; i < numPars; i++) {
+				System.arraycopy(products, i * numPars, h[i], 0, numPars);
+				h[i][i] += h[i][i] * lambda; // same LMA convention as getWJtJlambda()
+			}
+			final double [] b = new double [numPars];
+			System.arraycopy(products, numPars * numPars, b, 0, numPars);
+			wjtjlambda = new Matrix(h);
+			jty = new Matrix(b, b.length);
+		} else {
+			Matrix y_minus_fx_weighted = new Matrix(this.last_ymfx, this.last_ymfx.length);
+			wjtjlambda = new Matrix(getWJtJlambda(
+					lambda, // *10, // temporary
+					this.last_jt)); // double [][] jt)
+			jty = (new Matrix(this.last_jt)).times(y_minus_fx_weighted);
+		}
 		
 		if (debug_level>2) {
 			System.out.println("JtJ + lambda*diag(JtJ");
@@ -876,7 +909,6 @@ public class IntersceneLma {
 			jtjl_inv.print(18, 6);
 		}
 //last_jt has NaNs
-		Matrix jty = (new Matrix(this.last_jt)).times(y_minus_fx_weighted);
 		if (debug_level>2) {
 			System.out.println("Jt * (y-fx)");
 			jty.print(18, 6);