add a utest for operator_p_stack along the last second dimension
@@ -86,4 +86,55 @@ static bool test_op_p_stack(void)
UT_REGISTER_TEST(test_op_p_stack);
+static bool test_op_p_stack2(void)
+{
+ enum { N = 4 };
+ long dims[N] = { 8, 4, 1, 4 };
+ long dims2[N] = { 8, 4, 2, 4 };
+
+ long dims_no3[N] = { 8, 4, 1, 1 };
+
+ auto a = operator_p_scale(N, dims_no3);
+ auto a2 = operator_p_scale(N, dims_no3);
+
+ auto stack = operator_p_stack(2, 2, a, a2);
+
+ auto b = operator_p_stack(2, 2, a, a2);
+
+ long phases = dims[3];
+
+ for (int k = 0; k < (phases - 1); k++) {
+
+ auto tmp = operator_p_stack(3, 3, b, stack);
+ operator_p_free(b);
+ b = tmp;
+ }
+
+ operator_p_free(a);
+ operator_p_free(a2);
+ operator_p_free(stack);
+
+ complex float* in = md_alloc(N, dims2, CFL_SIZE);
+ complex float* out = md_alloc(N, dims2, CFL_SIZE);
+
+ md_zfill(N, dims2, in, 1.);
+ md_zfill(N, dims2, out, 100.);
+
+ operator_p_apply(b, 2., N, dims2, out, N, dims2, in);
+
+ operator_p_free(b);
+
+ md_zfill(N, dims2, in, 2.);
+
+ float err = md_znrmse(N, dims2, out, in);
+
+ md_free(in);
+ md_free(out);
+
+ return (err < UT_TOL);
+}
+
+UT_REGISTER_TEST(test_op_p_stack2);
+
+
