CODEX: Add deterministic LMA normal products kernel

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

jna/test_lma_products_jna.cu

+// test_lma_products_jna.cu - NVRTC/TpProc equivalence gate for roadmap rung 3.
+// It drives the same tp_proc_exec_lma_products entry point used by GpuQuadJna
+// and compares H/b bit-for-bit with a CPU oracle that enforces Java's
+// multiply/multiply/add expression order.
+//
+// Build after jna/build_lib.sh:
+//   cd tile_processor_gpu && /usr/local/cuda-12.8/bin/nvcc -std=c++17 -O2 \
+//     jna/test_lma_products_jna.cu -o tests_bin/test_lma_products_jna \
+//     -L jna -ltileproc -Xlinker -rpath -Xlinker $PWD/jna
+//
+// Created Jul 14, 2026 by Codex for Elphel.
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <vector>
+
+extern "C" {
+  void* tp_create_module(const char*, const char*);
+  const char* tp_last_error();
+  void* tp_proc_create(void*);
+  int   tp_proc_setup(void*, int,int,int,int,int,int);
+  int   tp_proc_exec_lma_products(void*, const double*, const double*, const double*, int,int,double*);
+  void  tp_proc_destroy(void*);
+  void  tp_destroy_module(void*);
+}
+
+static double mulAddJava(double sum, double a, double b)
+{
+  volatile double product = a * b;
+  volatile double next = sum + product;
+  return next;
+}
+
+int main(int argc, char** argv)
+{
+  const char* srcdir = (argc > 1) ? argv[1] : "src";
+  const char* devrt = (argc > 2) ? argv[2] :
+      "/usr/local/cuda-12.8/targets/x86_64-linux/lib/libcudadevrt.a";
+  void* module = tp_create_module(srcdir, devrt);
+  if(!module){ std::fprintf(stderr,"tp_create_module: %s\n",tp_last_error()); return EXIT_FAILURE; }
+  void* proc = tp_proc_create(module);
+  if(!proc || tp_proc_setup(proc,1,1,16,16,0,0)){
+    std::fprintf(stderr,"tp_proc setup: %s\n",tp_last_error());
+    if(proc) tp_proc_destroy(proc); tp_destroy_module(module); return EXIT_FAILURE;
+  }
+
+  const int np=3, nv=303; // current lean shape: 2*150 samples + 3 regularization rows
+  std::vector<double> weights(nv), jt((size_t)np*nv), ymfx(nv);
+  for(int k=0;k<nv;k++){
+    weights[k]=((k%11)==0)?0.0:(0.00037*(k+1)/(k+3.25));
+    ymfx[k]=((k&1)?-1.0:1.0)*(0.017+0.000031*k*k);
+    for(int i=0;i<np;i++){
+      const double sign=((k+2*i)&1)?-1.0:1.0;
+      jt[(size_t)i*nv+k]=sign*
+          (13.7*(i+1)+0.00091*(k+1)+0.0000013*((k*(i+3))%17));
+    }
+  }
+  std::vector<double> expected(np*np+np,0.0), got(expected.size(),0.0);
+  for(int i=0;i<np;i++){
+    for(int j=i;j<np;j++){
+      double sum=0.0;
+      for(int k=0;k<nv;k++){
+        volatile double wji=weights[k]*jt[(size_t)i*nv+k];
+        sum=mulAddJava(sum,wji,jt[(size_t)j*nv+k]);
+      }
+      expected[i*np+j]=expected[j*np+i]=sum;
+    }
+    double sum=0.0;
+    for(int k=0;k<nv;k++) sum=mulAddJava(sum,jt[(size_t)i*nv+k],ymfx[k]);
+    expected[np*np+i]=sum;
+  }
+
+  const int rc=tp_proc_exec_lma_products(proc,weights.data(),jt.data(),ymfx.data(),np,nv,got.data());
+  int bad=0;
+  if(rc){ std::fprintf(stderr,"tp_proc_exec_lma_products rc=%d: %s\n",rc,tp_last_error()); bad=1; }
+  else for(size_t i=0;i<got.size();i++) if(std::memcmp(&got[i],&expected[i],sizeof(double))!=0){
+    std::printf("[%zu] got %.17g expected %.17g -> FAIL\n",i,got[i],expected[i]); bad++;
+  }
+  std::printf("NVRTC lma_normal_products: %zu doubles, bit-exact -> %s\n",
+      got.size(),bad?"FAIL":"PASS");
+  tp_proc_destroy(proc); tp_destroy_module(module);
+  return bad?EXIT_FAILURE:EXIT_SUCCESS;
+}

src/tests/build_tests.sh

@@ -25,7 +25,8 @@ fi
 
 # Kernel + host sources shared by every test ("load all kernels"):
 COMMON="$SRC/TileProcessor.cu $SRC/dtt8x8.cu $SRC/geometry_correction.cu \
-        $SRC/tp_consolidate.cu $SRC/tp_peak.cu $SRC/tp_utils.cu $SRC/tp_files.cu tp_test_data.cu"
+        $SRC/tp_consolidate.cu $SRC/tp_peak.cu $SRC/tp_lma.cu \
+        $SRC/tp_utils.cu $SRC/tp_files.cu tp_test_data.cu"
 
 TESTS="${@:-$(ls test_*.cu | sed 's/\.cu$//')}"
 for t in $TESTS; do

src/tests/test_lma_products.cu

+/*
+ * test_lma_products.cu
+ *
+ * Standalone/cuda-gdb gate for lma_normal_products. The CPU oracle uses
+ * volatile intermediates to enforce the Java expression order and compares
+ * every output bit exactly (including a deliberately asymmetric sign/range
+ * mix and zero-weight samples).
+ *
+ * Usage: tests_bin/test_lma_products
+ *
+ * Created Jul 14, 2026 by Codex for Elphel (roadmap rung 3).
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <vector>
+
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
+#include "tp_lma.h"
+#include "tests/tp_gpu_buf.h"
+
+static double mulAddJava(double sum, double a, double b)
+{
+	volatile double product = a * b;
+	volatile double next = sum + product;
+	return next;
+}
+
+int main(int argc, char ** argv)
+{
+	findCudaDevice(argc, (const char **) argv);
+	const int np = 3;
+	const int nv = 303; // current lean shape: 2*150 samples + 3 regularization rows
+	std::vector<double> weights(nv), jt((size_t) np * nv), ymfx(nv);
+	for (int k = 0; k < nv; k++) {
+		weights[k] = ((k % 11) == 0) ? 0.0 : (0.00037 * (k + 1) / (k + 3.25));
+		ymfx[k] = ((k & 1) ? -1.0 : 1.0) * (0.017 + 0.000031 * k * k);
+		for (int i = 0; i < np; i++) {
+			const double sign = ((k + 2 * i) & 1) ? -1.0 : 1.0;
+			jt[(size_t) i * nv + k] = sign *
+					(13.7 * (i + 1) + 0.00091 * (k + 1) + 0.0000013 * ((k * (i + 3)) % 17));
+		}
+	}
+
+	std::vector<double> expected(np * np + np, 0.0);
+	for (int i = 0; i < np; i++) {
+		for (int j = i; j < np; j++) {
+			double sum = 0.0;
+			for (int k = 0; k < nv; k++) {
+				volatile double wji = weights[k] * jt[(size_t) i * nv + k];
+				sum = mulAddJava(sum, wji, jt[(size_t) j * nv + k]);
+			}
+			expected[i * np + j] = expected[j * np + i] = sum;
+		}
+		double sum = 0.0;
+		for (int k = 0; k < nv; k++) {
+			sum = mulAddJava(sum, jt[(size_t) i * nv + k], ymfx[k]);
+		}
+		expected[np * np + i] = sum;
+	}
+
+	GpuBuf<double> gpu_w(nv), gpu_jt((size_t) np * nv), gpu_y(nv), gpu_out(expected.size());
+	gpu_w.upload(weights.data());
+	gpu_jt.upload(jt.data());
+	gpu_y.upload(ymfx.data());
+	lma_normal_products<<<1, 1>>>(np, nv, gpu_w.dev(), gpu_jt.dev(), gpu_y.dev(), gpu_out.dev());
+	checkCudaErrors(cudaGetLastError());
+	checkCudaErrors(cudaDeviceSynchronize());
+	const std::vector<double> got = gpu_out.download();
+
+	int bad = 0;
+	for (size_t i = 0; i < got.size(); i++) {
+		if (memcmp(&got[i], &expected[i], sizeof(double)) != 0) {
+			printf("[%zu] got %.17g expected %.17g -> FAIL\n", i, got[i], expected[i]);
+			bad++;
+		}
+	}
+	printf("lma_normal_products: %zu doubles, bit-exact -> %s\n",
+			got.size(), bad ? "FAIL" : "PASS");
+	return bad ? EXIT_FAILURE : EXIT_SUCCESS;
+}

src/tp_lma.cu

+/*
+ * tp_lma.cu
+ *
+ * Deterministic double-precision normal-equation products for the lean pose
+ * LMA. One thread intentionally follows IntersceneLma's Java summation order:
+ * the arrays are tiny (currently 3 parameters x roughly 303 values), while a
+ * parallel reduction would change rounding and undermine the tol-0 validation
+ * tier. Explicit __dmul_rn/__dadd_rn prevent FMA contraction from merging the
+ * Java oracle's two multiplies and add.
+ *
+ * Created Jul 14, 2026 by Codex for Elphel (roadmap rung 3).
+ */
+
+#ifndef JCUDA
+#include <cuda_runtime.h>
+#include "tp_lma.h"
+#endif
+
+extern "C" __global__ void lma_normal_products(
+		int            num_params,
+		int            num_values,
+		const double * weights,
+		const double * jt,
+		const double * ymfx_weighted,
+		double       * products)
+{
+	if ((blockIdx.x != 0) || (threadIdx.x != 0)) return;
+
+	for (int i = 0; i < num_params; i++) {
+		const double * ji = jt + (size_t) i * num_values;
+		for (int j = i; j < num_params; j++) {
+			const double * jj = jt + (size_t) j * num_values;
+			double sum = 0.0;
+			for (int k = 0; k < num_values; k++) {
+				const double wji = __dmul_rn(weights[k], ji[k]);
+				const double term = __dmul_rn(wji, jj[k]);
+				sum = __dadd_rn(sum, term);
+			}
+			products[i * num_params + j] = sum;
+			products[j * num_params + i] = sum;
+		}
+
+		double sum = 0.0;
+		for (int k = 0; k < num_values; k++) {
+			const double term = __dmul_rn(ji[k], ymfx_weighted[k]);
+			sum = __dadd_rn(sum, term);
+		}
+		products[num_params * num_params + i] = sum;
+	}
+}

src/tp_lma.h

+/*
+ * tp_lma.h
+ *
+ * Small normal-equation product kernel for the pose-chain LMA. Host policy
+ * chooses the active parameters/tiles and prepares the eigen-transformed
+ * Jacobian, weights and weighted residual; the GPU returns only
+ *
+ *   H = J^T W J    and    b = J^T W (y - f(x)).
+ *
+ * The kernel contract is parameterized by num_params/num_values. The current
+ * CUAS lean path is 3 parameters and about 2*150 values, but the kernel does
+ * not encode those policy constants.
+ *
+ * Created Jul 14, 2026 by Codex for Elphel (roadmap rung 3).
+ */
+
+#ifndef SRC_TP_LMA_H_
+#define SRC_TP_LMA_H_
+
+#ifndef JCUDA
+#include <cuda_runtime.h>
+#endif
+
+extern "C" __global__ void lma_normal_products(
+		int            num_params,
+		int            num_values,
+		const double * weights,
+		const double * jt,             // parameter-major [num_params][num_values]
+		const double * ymfx_weighted,  // W * (y - f), already eigen-transformed
+		double       * products);      // row-major H [p*p], then b [p]
+
+#endif /* SRC_TP_LMA_H_ */