CODEX: Add float CUDA pose LMA candidate

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

jna/test_pose_lma_step_jna.cu

+// test_pose_lma_step_jna.cu - NVRTC/TpProc equivalence gate for rung 3-A4d.
+// Drives the same H2D/kernel/D2H entry point used by GpuQuadJna and compares
+// every float against a CPU oracle with Java-compatible operation ordering.
+//
+// Build after jna/build_lib.sh:
+//   cd tile_processor_gpu && /usr/local/cuda-12.8/bin/nvcc -std=c++17 -O2 \
+//     jna/test_pose_lma_step_jna.cu -o tests_bin/test_pose_lma_step_jna \
+//     -L jna -ltileproc -Xlinker -rpath -Xlinker $PWD/jna
+// Run: tests_bin/test_pose_lma_step_jna [srcdir] [libcudadevrt.a]
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <vector>
+
+enum { NUM_PARAMS = 3, PRODUCTS = 12, RESULT_FLOATS = 19 };
+
+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_pose_lma_step(void*, float, const float*, const float*,
+      const float*, const float*, int, float*);
+  void  tp_proc_destroy(void*);
+  void  tp_destroy_module(void*);
+}
+
+static float fmul(float a, float b) { volatile float v = a * b; return v; }
+static float fadd(float a, float b) { volatile float v = a + b; return v; }
+static float fsub(float a, float b) { volatile float v = a - b; return v; }
+static float fdiv(float a, float b) { volatile float v = a / b; return v; }
+
+static void cpuPoseLmaStep(int nv, float lambda, const float* weights,
+    const float* jt, const float* ymfx, const float* vector, float* result)
+{
+  result[0]=0.0f;
+  for(int i=1;i<RESULT_FLOATS;i++) result[i]=NAN;
+  for(int i=0;i<NUM_PARAMS;i++) result[16+i]=vector[i];
+  float products[PRODUCTS]={0.0f};
+  for(int row=0;row<NUM_PARAMS;row++){
+    const float* jr=jt+(size_t)row*nv;
+    for(int column=row;column<NUM_PARAMS;column++){
+      const float* jc=jt+(size_t)column*nv;
+      float sum=0.0f;
+      for(int value=0;value<nv;value++){
+        const float weighted=fmul(weights[value],jr[value]);
+        sum=fadd(sum,fmul(weighted,jc[value]));
+      }
+      products[row*NUM_PARAMS+column]=products[column*NUM_PARAMS+row]=sum;
+    }
+    float sum=0.0f;
+    for(int value=0;value<nv;value++) sum=fadd(sum,fmul(jr[value],ymfx[value]));
+    products[9+row]=sum;
+  }
+  for(int i=0;i<PRODUCTS;i++) result[1+i]=products[i];
+  float a[9]; for(int i=0;i<9;i++) a[i]=products[i];
+  for(int i=0;i<NUM_PARAMS;i++){
+    const int diagonal=i*NUM_PARAMS+i;
+    a[diagonal]=fadd(a[diagonal],fmul(a[diagonal],lambda));
+  }
+  const float c00=fsub(fmul(a[4],a[8]),fmul(a[5],a[7]));
+  const float c01=fsub(fmul(a[2],a[7]),fmul(a[1],a[8]));
+  const float c02=fsub(fmul(a[1],a[5]),fmul(a[2],a[4]));
+  const float c10=fsub(fmul(a[5],a[6]),fmul(a[3],a[8]));
+  const float c11=fsub(fmul(a[0],a[8]),fmul(a[2],a[6]));
+  const float c12=fsub(fmul(a[2],a[3]),fmul(a[0],a[5]));
+  const float c20=fsub(fmul(a[3],a[7]),fmul(a[4],a[6]));
+  const float c21=fsub(fmul(a[1],a[6]),fmul(a[0],a[7]));
+  const float c22=fsub(fmul(a[0],a[4]),fmul(a[1],a[3]));
+  float det=fadd(fmul(a[0],c00),fmul(a[1],c10));
+  det=fadd(det,fmul(a[2],c20));
+  if(!std::isfinite(det)||(det==0.0f)) return;
+  const float scale=fdiv(1.0f,det);
+  const float inverse[9]={
+      fmul(c00,scale),fmul(c01,scale),fmul(c02,scale),
+      fmul(c10,scale),fmul(c11,scale),fmul(c12,scale),
+      fmul(c20,scale),fmul(c21,scale),fmul(c22,scale)};
+  for(int row=0;row<NUM_PARAMS;row++){
+    float delta=fadd(fmul(inverse[3*row],products[9]),fmul(inverse[3*row+1],products[10]));
+    delta=fadd(delta,fmul(inverse[3*row+2],products[11]));
+    result[13+row]=delta; result[16+row]=fadd(vector[row],delta);
+  }
+  for(int row=0;row<NUM_PARAMS;row++)
+    if(!std::isfinite(result[13+row])||!std::isfinite(result[16+row])) return;
+  result[0]=1.0f;
+}
+
+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 nv=303; const float lambda=0.001f;
+  std::vector<float> weights(nv),jt((size_t)NUM_PARAMS*nv),ymfx(nv);
+  const float vector[NUM_PARAMS]={0.0023f,-0.0017f,0.00091f};
+  for(int k=0;k<nv;k++){
+    weights[k]=((k%11)==0)?0.0f:(0.00037f*(k+1)/(k+3.25f));
+    ymfx[k]=((k&1)?-1.0f:1.0f)*(0.017f+0.000031f*k*k);
+    for(int i=0;i<NUM_PARAMS;i++){
+      const float sign=((k+3*i)&1)?-1.0f:1.0f;
+      jt[(size_t)i*nv+k]=sign*(13.7f*(i+1)+0.00091f*(k+1)+
+          0.0000013f*((k*(i+3))%17));
+    }
+  }
+  float expected[RESULT_FLOATS],got[RESULT_FLOATS];
+  cpuPoseLmaStep(nv,lambda,weights.data(),jt.data(),ymfx.data(),vector,expected);
+  const int rc=tp_proc_exec_pose_lma_step(proc,lambda,weights.data(),jt.data(),
+      ymfx.data(),vector,nv,got);
+  int bad=0;
+  if(rc){ std::fprintf(stderr,"tp_proc_exec_pose_lma_step rc=%d: %s\n",rc,tp_last_error()); bad=1; }
+  else for(int i=0;i<RESULT_FLOATS;i++) if(std::memcmp(&got[i],&expected[i],sizeof(float))!=0){
+    std::printf("[%d] got %.9g expected %.9g -> FAIL\n",i,got[i],expected[i]); bad++;
+  }
+  std::printf("NVRTC pose_lma_step: %d floats, bit-exact -> %s\n",
+      RESULT_FLOATS,bad?"FAIL":"PASS");
+  tp_proc_destroy(proc); tp_destroy_module(module);
+  return bad?EXIT_FAILURE:EXIT_SUCCESS;
+}

jna/tp_jna.cpp

@@ -54,6 +54,7 @@ static const char* KERNELS[] = {
   "clt_average_sensors","index_consolidate","consolidate_oob","clt_average_sensors_list", // By Claude on 07/12/2026
   "corr2D_peak_eig","corr2D_peak_eig_f", // double + float instantiations (precision two-step) // By Claude on 07/13-14/2026
   "pose_fx_jacobian", // one float thread per selected tile // By Codex on 07/15/2026
+  "pose_lma_step", // fixed 3-angle float products + direct solve in one thread // By Codex on 07/15/2026
   "lma_normal_products"}; // deterministic double normal-equation products // By Codex on 07/14/2026
 static const int N_KERNELS = sizeof(KERNELS)/sizeof(KERNELS[0]);
 
@@ -654,6 +655,11 @@ struct TpProc {
   float *gpu_pose_radial[2], *gpu_pose_rbr[2], *gpu_pose_ers[2];
   float *gpu_pose_vectors, *gpu_pose_centers, *gpu_pose_fx, *gpu_pose_jt;
   unsigned char *gpu_pose_selection, *gpu_pose_valid;
+  // One-thread float runLma candidate validation rung. Host-prepared J/W/residual
+  // arrays intentionally cross H2D and the compact result crosses D2H.
+  int pose_lma_value_capacity;
+  float *gpu_pose_lma_weights, *gpu_pose_lma_jt, *gpu_pose_lma_ymfx;
+  float *gpu_pose_lma_vector, *gpu_pose_lma_result;
 };
 
 extern "C" {
@@ -695,6 +701,9 @@ TpProc* tp_proc_create(TpModule* m){
   p->gpu_pose_ers[0]=p->gpu_pose_ers[1]=nullptr;
   p->gpu_pose_vectors=p->gpu_pose_centers=p->gpu_pose_fx=p->gpu_pose_jt=nullptr;
   p->gpu_pose_selection=p->gpu_pose_valid=nullptr; // By Codex on 07/15/2026
+  p->pose_lma_value_capacity=0;
+  p->gpu_pose_lma_weights=p->gpu_pose_lma_jt=p->gpu_pose_lma_ymfx=nullptr;
+  p->gpu_pose_lma_vector=p->gpu_pose_lma_result=nullptr; // By Codex on 07/15/2026
   return p;
 }
 
@@ -1295,6 +1304,50 @@ int tp_proc_exec_pose_fx_jacobian(TpProc* p,
   return 0;
 }
 
+// Primitive-float fixed three-angle candidate calculation. This validation
+// boundary deliberately uploads transformed J/W/residual state on every call;
+// the final resident chain will produce and consume those arrays on-device.
+int tp_proc_exec_pose_lma_step(TpProc* p,
+                               float lambda,
+                               const float* weights,
+                               const float* jt,
+                               const float* ymfx_weighted,
+                               const float* vector,
+                               int num_values,
+                               float* result){
+  if(!p||!weights||!jt||!ymfx_weighted||!vector||!result||(num_values<=0)){
+    seterr("exec_pose_lma_step: bad args p=%p nv=%d",(void*)p,num_values); return -1; }
+  cuCtxSetCurrent(p->mod->ctx);
+  CUfunction f=getfun(p->mod,"pose_lma_step");
+  if(!f){ seterr("pose_lma_step missing"); return -2; }
+  if(num_values>p->pose_lma_value_capacity){
+    cudaFree(p->gpu_pose_lma_weights); cudaFree(p->gpu_pose_lma_jt);
+    cudaFree(p->gpu_pose_lma_ymfx); cudaFree(p->gpu_pose_lma_vector);
+    cudaFree(p->gpu_pose_lma_result);
+    p->gpu_pose_lma_weights=p->gpu_pose_lma_jt=p->gpu_pose_lma_ymfx=nullptr;
+    p->gpu_pose_lma_vector=p->gpu_pose_lma_result=nullptr;
+    p->pose_lma_value_capacity=0;
+    if(cudaMalloc((void**)&p->gpu_pose_lma_weights,(size_t)num_values*sizeof(float))!=cudaSuccess ||
+       cudaMalloc((void**)&p->gpu_pose_lma_jt,(size_t)POSE_NUM_PARAMS*num_values*sizeof(float))!=cudaSuccess ||
+       cudaMalloc((void**)&p->gpu_pose_lma_ymfx,(size_t)num_values*sizeof(float))!=cudaSuccess ||
+       cudaMalloc((void**)&p->gpu_pose_lma_vector,POSE_NUM_PARAMS*sizeof(float))!=cudaSuccess ||
+       cudaMalloc((void**)&p->gpu_pose_lma_result,POSE_LMA_RESULT*sizeof(float))!=cudaSuccess){
+      seterr("exec_pose_lma_step: scratch alloc failed nv=%d",num_values); return -3; }
+    p->pose_lma_value_capacity=num_values;
+  }
+  if(cudaMemcpy(p->gpu_pose_lma_weights,weights,(size_t)num_values*sizeof(float),cudaMemcpyHostToDevice)!=cudaSuccess ||
+     cudaMemcpy(p->gpu_pose_lma_jt,jt,(size_t)POSE_NUM_PARAMS*num_values*sizeof(float),cudaMemcpyHostToDevice)!=cudaSuccess ||
+     cudaMemcpy(p->gpu_pose_lma_ymfx,ymfx_weighted,(size_t)num_values*sizeof(float),cudaMemcpyHostToDevice)!=cudaSuccess ||
+     cudaMemcpy(p->gpu_pose_lma_vector,vector,POSE_NUM_PARAMS*sizeof(float),cudaMemcpyHostToDevice)!=cudaSuccess){
+    seterr("exec_pose_lma_step: HtoD failed"); return -4; }
+  void* a[]={ &num_values,&lambda,&p->gpu_pose_lma_weights,&p->gpu_pose_lma_jt,
+              &p->gpu_pose_lma_ymfx,&p->gpu_pose_lma_vector,&p->gpu_pose_lma_result };
+  if(launch1(f,1,1,1,1,1,1,a,"pose_lma_step")) return -5;
+  if(cudaMemcpy(result,p->gpu_pose_lma_result,POSE_LMA_RESULT*sizeof(float),cudaMemcpyDeviceToHost)!=cudaSuccess){
+    seterr("exec_pose_lma_step: DtoH failed"); return -6; }
+  return 0;
+}
+
 // Deterministic double normal-equation products for the lean pose LMA:
 // products = row-major H[num_params][num_params], followed by b[num_params].
 // The input arrays are already eigen-transformed/weighted by IntersceneLma;
@@ -1549,6 +1602,8 @@ void tp_proc_destroy(TpProc* p){
   cudaFree(p->gpu_pose_rbr[0]); cudaFree(p->gpu_pose_rbr[1]); cudaFree(p->gpu_pose_ers[0]); cudaFree(p->gpu_pose_ers[1]);
   cudaFree(p->gpu_pose_vectors); cudaFree(p->gpu_pose_centers); cudaFree(p->gpu_pose_selection);
   cudaFree(p->gpu_pose_fx); cudaFree(p->gpu_pose_jt); cudaFree(p->gpu_pose_valid); // By Codex on 07/15/2026
+  cudaFree(p->gpu_pose_lma_weights); cudaFree(p->gpu_pose_lma_jt); cudaFree(p->gpu_pose_lma_ymfx);
+  cudaFree(p->gpu_pose_lma_vector); cudaFree(p->gpu_pose_lma_result); // By Codex on 07/15/2026
   delete p;
 }
 

src/tests/test_pose_lma_step.cu

+/*
+ * test_pose_lma_step.cu
+ *
+ * Standalone bit-exact gate for the fixed three-angle float LMA candidate.
+ * The CPU oracle forces a float rounding after every elementary operation so
+ * it has the same expression order as Java and the CUDA __f*_rn operations.
+ *
+ * Usage: tests_bin/test_pose_lma_step
+ *
+ * Created Jul 15, 2026 by Codex for Elphel (roadmap rung 3-A4d).
+ */
+
+#include <cmath>
+#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 float fmul(float a, float b) { volatile float v = a * b; return v; }
+static float fadd(float a, float b) { volatile float v = a + b; return v; }
+static float fsub(float a, float b) { volatile float v = a - b; return v; }
+static float fdiv(float a, float b) { volatile float v = a / b; return v; }
+
+static void cpuPoseLmaStep(
+		int numValues,
+		float lambda,
+		const float * weights,
+		const float * jt,
+		const float * ymfx,
+		const float * vector,
+		float * result)
+{
+	result[0] = 0.0f;
+	for (int i = 1; i < POSE_LMA_RESULT; i++) result[i] = NAN;
+	for (int i = 0; i < POSE_NUM_PARAMS; i++) result[16 + i] = vector[i];
+	float products[POSE_LMA_PRODUCTS] = {0.0f};
+	for (int row = 0; row < POSE_NUM_PARAMS; row++) {
+		const float * jr = jt + (size_t) row * numValues;
+		for (int column = row; column < POSE_NUM_PARAMS; column++) {
+			const float * jc = jt + (size_t) column * numValues;
+			float sum = 0.0f;
+			for (int value = 0; value < numValues; value++) {
+				const float weighted = fmul(weights[value], jr[value]);
+				sum = fadd(sum, fmul(weighted, jc[value]));
+			}
+			products[row * POSE_NUM_PARAMS + column] = sum;
+			products[column * POSE_NUM_PARAMS + row] = sum;
+		}
+		float sum = 0.0f;
+		for (int value = 0; value < numValues; value++) {
+			sum = fadd(sum, fmul(jr[value], ymfx[value]));
+		}
+		products[9 + row] = sum;
+	}
+	for (int i = 0; i < POSE_LMA_PRODUCTS; i++) result[1 + i] = products[i];
+	float a[9];
+	for (int i = 0; i < 9; i++) a[i] = products[i];
+	for (int i = 0; i < POSE_NUM_PARAMS; i++) {
+		const int diagonal = i * POSE_NUM_PARAMS + i;
+		a[diagonal] = fadd(a[diagonal], fmul(a[diagonal], lambda));
+	}
+	const float c00 = fsub(fmul(a[4],a[8]), fmul(a[5],a[7]));
+	const float c01 = fsub(fmul(a[2],a[7]), fmul(a[1],a[8]));
+	const float c02 = fsub(fmul(a[1],a[5]), fmul(a[2],a[4]));
+	const float c10 = fsub(fmul(a[5],a[6]), fmul(a[3],a[8]));
+	const float c11 = fsub(fmul(a[0],a[8]), fmul(a[2],a[6]));
+	const float c12 = fsub(fmul(a[2],a[3]), fmul(a[0],a[5]));
+	const float c20 = fsub(fmul(a[3],a[7]), fmul(a[4],a[6]));
+	const float c21 = fsub(fmul(a[1],a[6]), fmul(a[0],a[7]));
+	const float c22 = fsub(fmul(a[0],a[4]), fmul(a[1],a[3]));
+	float determinant = fadd(fmul(a[0],c00), fmul(a[1],c10));
+	determinant = fadd(determinant, fmul(a[2],c20));
+	if (!std::isfinite(determinant) || (determinant == 0.0f)) return;
+	const float scale = fdiv(1.0f, determinant);
+	const float inverse[9] = {
+			fmul(c00,scale),fmul(c01,scale),fmul(c02,scale),
+			fmul(c10,scale),fmul(c11,scale),fmul(c12,scale),
+			fmul(c20,scale),fmul(c21,scale),fmul(c22,scale)};
+	for (int row = 0; row < POSE_NUM_PARAMS; row++) {
+		float delta = fadd(fmul(inverse[3 * row], products[9]),
+				fmul(inverse[3 * row + 1], products[10]));
+		delta = fadd(delta, fmul(inverse[3 * row + 2], products[11]));
+		result[13 + row] = delta;
+		result[16 + row] = fadd(vector[row], delta);
+	}
+	for (int row = 0; row < POSE_NUM_PARAMS; row++)
+		if (!std::isfinite(result[13 + row]) || !std::isfinite(result[16 + row])) return;
+	result[0] = 1.0f;
+}
+
+int main(int argc, char ** argv)
+{
+	findCudaDevice(argc, (const char **) argv);
+	const int nv = 303;
+	const float lambda = 0.001f;
+	std::vector<float> weights(nv), jt((size_t) POSE_NUM_PARAMS * nv), ymfx(nv);
+	const float vector[POSE_NUM_PARAMS] = {0.0023f, -0.0017f, 0.00091f};
+	for (int k = 0; k < nv; k++) {
+		weights[k] = ((k % 11) == 0) ? 0.0f : (0.00037f * (k + 1) / (k + 3.25f));
+		ymfx[k] = ((k & 1) ? -1.0f : 1.0f) * (0.017f + 0.000031f * k * k);
+		for (int i = 0; i < POSE_NUM_PARAMS; i++) {
+			const float sign = ((k + 3 * i) & 1) ? -1.0f : 1.0f;
+			jt[(size_t) i * nv + k] = sign *
+					(13.7f * (i + 1) + 0.00091f * (k + 1) +
+					0.0000013f * ((k * (i + 3)) % 17));
+		}
+	}
+	float expected[POSE_LMA_RESULT];
+	cpuPoseLmaStep(nv, lambda, weights.data(), jt.data(), ymfx.data(), vector, expected);
+	GpuBuf<float> gpuWeights(nv), gpuJt((size_t) POSE_NUM_PARAMS * nv), gpuYmfx(nv);
+	GpuBuf<float> gpuVector(POSE_NUM_PARAMS), gpuResult(POSE_LMA_RESULT);
+	gpuWeights.upload(weights.data());
+	gpuJt.upload(jt.data());
+	gpuYmfx.upload(ymfx.data());
+	gpuVector.upload(vector);
+	pose_lma_step<<<1, 1>>>(nv, lambda, gpuWeights.dev(), gpuJt.dev(), gpuYmfx.dev(),
+			gpuVector.dev(), gpuResult.dev());
+	checkCudaErrors(cudaGetLastError());
+	checkCudaErrors(cudaDeviceSynchronize());
+	const std::vector<float> got = gpuResult.download();
+	int bad = 0;
+	for (int i = 0; i < POSE_LMA_RESULT; i++) {
+		if (std::memcmp(&got[i], &expected[i], sizeof(float)) != 0) {
+			std::printf("[%d] got %.9g expected %.9g -> FAIL\n", i, got[i], expected[i]);
+			bad++;
+		}
+	}
+	std::printf("pose_lma_step: %d floats, bit-exact -> %s\n",
+			POSE_LMA_RESULT, bad ? "FAIL" : "PASS");
+	return bad ? EXIT_FAILURE : EXIT_SUCCESS;
+}

src/tp_lma.cu

@@ -473,6 +473,102 @@ extern "C" __global__ void pose_fx_jacobian(
 	valid[tile] = 1;
 }
 
+static __device__ __forceinline__ float pose_lma_mul(float a, float b)
+{
+	return __fmul_rn(a, b);
+}
+
+static __device__ __forceinline__ float pose_lma_add(float a, float b)
+{
+	return __fadd_rn(a, b);
+}
+
+static __device__ __forceinline__ float pose_lma_sub(float a, float b)
+{
+	return __fsub_rn(a, b);
+}
+
+extern "C" __global__ void pose_lma_step(
+		int           num_values,
+		float         lambda,
+		const float * weights,
+		const float * jt,
+		const float * ymfx_weighted,
+		const float * vector,
+		float       * result)
+{
+	if ((blockIdx.x != 0) || (threadIdx.x != 0)) return;
+	result[0] = 0.0f;
+	#pragma unroll
+	for (int i = 1; i < 1 + POSE_LMA_PRODUCTS + POSE_NUM_PARAMS; i++) result[i] = nanf("");
+	#pragma unroll
+	for (int i = 0; i < POSE_NUM_PARAMS; i++) result[16 + i] = vector[i];
+
+	float products[POSE_LMA_PRODUCTS] = {0.0f};
+	#pragma unroll
+	for (int row = 0; row < POSE_NUM_PARAMS; row++) {
+		const float * jr = jt + (size_t) row * num_values;
+		#pragma unroll
+		for (int column = row; column < POSE_NUM_PARAMS; column++) {
+			const float * jc = jt + (size_t) column * num_values;
+			float sum = 0.0f;
+			for (int value = 0; value < num_values; value++) {
+				const float weighted = pose_lma_mul(weights[value], jr[value]);
+				sum = pose_lma_add(sum, pose_lma_mul(weighted, jc[value]));
+			}
+			products[row * POSE_NUM_PARAMS + column] = sum;
+			products[column * POSE_NUM_PARAMS + row] = sum;
+		}
+		float sum = 0.0f;
+		for (int value = 0; value < num_values; value++) {
+			sum = pose_lma_add(sum, pose_lma_mul(jr[value], ymfx_weighted[value]));
+		}
+		products[POSE_NUM_PARAMS * POSE_NUM_PARAMS + row] = sum;
+	}
+	#pragma unroll
+	for (int i = 0; i < POSE_LMA_PRODUCTS; i++) result[1 + i] = products[i];
+
+	float a[9];
+	#pragma unroll
+	for (int i = 0; i < 9; i++) a[i] = products[i];
+	#pragma unroll
+	for (int i = 0; i < POSE_NUM_PARAMS; i++) {
+		const int diagonal = i * POSE_NUM_PARAMS + i;
+		a[diagonal] = pose_lma_add(a[diagonal], pose_lma_mul(a[diagonal], lambda));
+	}
+	const float c00 = pose_lma_sub(pose_lma_mul(a[4],a[8]), pose_lma_mul(a[5],a[7]));
+	const float c01 = pose_lma_sub(pose_lma_mul(a[2],a[7]), pose_lma_mul(a[1],a[8]));
+	const float c02 = pose_lma_sub(pose_lma_mul(a[1],a[5]), pose_lma_mul(a[2],a[4]));
+	const float c10 = pose_lma_sub(pose_lma_mul(a[5],a[6]), pose_lma_mul(a[3],a[8]));
+	const float c11 = pose_lma_sub(pose_lma_mul(a[0],a[8]), pose_lma_mul(a[2],a[6]));
+	const float c12 = pose_lma_sub(pose_lma_mul(a[2],a[3]), pose_lma_mul(a[0],a[5]));
+	const float c20 = pose_lma_sub(pose_lma_mul(a[3],a[7]), pose_lma_mul(a[4],a[6]));
+	const float c21 = pose_lma_sub(pose_lma_mul(a[1],a[6]), pose_lma_mul(a[0],a[7]));
+	const float c22 = pose_lma_sub(pose_lma_mul(a[0],a[4]), pose_lma_mul(a[1],a[3]));
+	float determinant = pose_lma_add(pose_lma_mul(a[0],c00), pose_lma_mul(a[1],c10));
+	determinant = pose_lma_add(determinant, pose_lma_mul(a[2],c20));
+	if (!isfinite(determinant) || (determinant == 0.0f)) return;
+	const float scale = __fdiv_rn(1.0f, determinant);
+	const float inverse[9] = {
+		pose_lma_mul(c00,scale),pose_lma_mul(c01,scale),pose_lma_mul(c02,scale),
+		pose_lma_mul(c10,scale),pose_lma_mul(c11,scale),pose_lma_mul(c12,scale),
+		pose_lma_mul(c20,scale),pose_lma_mul(c21,scale),pose_lma_mul(c22,scale)};
+	const float * b = products + 9;
+	float delta[3];
+	#pragma unroll
+	for (int row = 0; row < POSE_NUM_PARAMS; row++) {
+		float value = pose_lma_add(
+				pose_lma_mul(inverse[3 * row], b[0]),
+				pose_lma_mul(inverse[3 * row + 1], b[1]));
+		value = pose_lma_add(value, pose_lma_mul(inverse[3 * row + 2], b[2]));
+		delta[row] = value;
+		result[13 + row] = value;
+		result[16 + row] = pose_lma_add(vector[row], value);
+	}
+	if (!pose_finite3(delta) || !pose_finite3(result + 16)) return;
+	result[0] = 1.0f;
+}
+
 extern "C" __global__ void lma_normal_products(
 		int            num_params,
 		int            num_values,

src/tp_lma.h

@@ -27,6 +27,8 @@
 #define POSE_NUM_COMPONENTS   2
 #define POSE_RADIAL_FLOATS    7
 #define POSE_ERS_LINE_STRIDE 14
+#define POSE_LMA_PRODUCTS    12
+#define POSE_LMA_RESULT      19
 
 typedef struct {
 	int   width;
@@ -57,6 +59,17 @@ extern "C" __global__ void pose_fx_jacobian(
 		float                * jt,             // parameter-major [3][2*num_tiles], nullable
 		unsigned char        * valid);
 
+// Fixed three-angle candidate calculation in one float thread. Output layout:
+// {valid, raw H[9], raw b[3], delta[3], candidate[3]}.
+extern "C" __global__ void pose_lma_step(
+		int           num_values,
+		float         lambda,
+		const float * weights,
+		const float * jt,             // parameter-major [3][num_values]
+		const float * ymfx_weighted,
+		const float * vector,         // current three-angle parameter vector
+		float       * result);        // [POSE_LMA_RESULT]
+
 extern "C" __global__ void lma_normal_products(
 		int            num_params,
 		int            num_values,