From af8b5739f1b60be08bf7de2ba6b2c02c1247f0ee Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Sat, 1 Aug 2026 02:12:29 -0700 Subject: [PATCH 1/3] fix: seed scm[0] before symmetry-mode sum overwrites b[i] solves() in src/matrix_algebra.cpp transforms the right hand side into symmetry modes established by a GX card, solves each mode, then transforms back. Both transforms seed each image row from the row's untransformed value, held in scm[0], but the accumulation loops folded the k=0 term directly into the running sum and began at k=1, so scm[0] was never written while the mode loops below continued to read it. scm[0] is allocated via resize(), which does not zero, so wire-only structures with a GX card produced subnormal doubles near 1e-311 in currents and INF in impedance rather than a fault. - store scm[0] = b[i+column_offset] before the accumulation overwrites b[i], in both the forward and inverse mode transforms, matching nec2c.c:10523 and xnec2c/src/matrix.c:1490 Signed-off-by: Eric Wheeler --- src/matrix_algebra.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/matrix_algebra.cpp b/src/matrix_algebra.cpp index 817b760..6f9ff27 100644 --- a/src/matrix_algebra.cpp +++ b/src/matrix_algebra.cpp @@ -398,7 +398,11 @@ void solves(complex_array& a, int_array& ip, complex_array& b, int64_t neq, /* transform matrix eq. rhs vector according to symmetry modes */ for (int64_t i = 0; i < npeq; i++ ) { - nec_complex sum_normal(b[i+column_offset]); + /* Each image row below is seeded from this row's untransformed + value, so hold it in scm[0] before the normalized sum + overwrites b[i]. */ + scm[0] = b[i+column_offset]; + nec_complex sum_normal(scm[0]); for (int64_t k = 1; k < nop; k++ ) { int64_t ia= i+ k* npeq; scm[k]= b[ia+column_offset]; @@ -441,7 +445,11 @@ void solves(complex_array& a, int_array& ip, complex_array& b, int64_t neq, for (int64_t ic = 0; ic < nrh; ic++ ) { int64_t column_offset = ic*neq; for (int64_t i = 0; i < npeq; i++ ) { - nec_complex sum_normal(b[i+column_offset]); + /* Each image row below is seeded from this row's untransformed + value, so hold it in scm[0] before the summed value overwrites + b[i]. */ + scm[0] = b[i+column_offset]; + nec_complex sum_normal(scm[0]); for (int64_t k = 1; k < nop; k++ ) { int64_t ia= i+ k* npeq; scm[k]= b[ia+column_offset]; From 33f93d89022cf32fdeb64521d5828ac43aa388fc Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Sat, 1 Aug 2026 02:19:25 -0700 Subject: [PATCH 2/3] fix: correct plane-symmetry matrix fill guard in fblock() nec_context::fblock() builds symmetry_array, the matrix of modal coefficients used to expand a GX-reduced solve back onto the image sections. The Fortran this routine transcribes tests the section count and falls through to the fill loop only for a legal nop of 2, 4, or 8, otherwise stopping. The C++ inverted that sense: it returned early for the legal counts and reached the fill loop only for an illegal one, so every symmetric structure left symmetry_array populated at element 0 alone, the rest uninitialized. The pass count was also computed as nop / 2 rather than log2(nop), which agrees with the correct value at 2 and 4 but diverges at 8. - restore the Fortran's guard sense, stopping on an illegal nop instead of returning on a legal one, so the fill loop always executes - derive the pass count by doubling until it reaches nop, matching nec2c.c's derivation, instead of nop / 2 Combined with the prior scm[0] seeding fix, this closes the uninitialized symmetry_array read that produced subnormal currents and infinite impedance on every GX-symmetric structure. Signed-off-by: Eric Wheeler --- src/nec_context.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/nec_context.cpp b/src/nec_context.cpp index 6edf200..df5faf6 100644 --- a/src/nec_context.cpp +++ b/src/nec_context.cpp @@ -6439,14 +6439,17 @@ void nec_context::fblock( int nrow, int ncol, int64_t imax, int ipsym ) { int kk=1; symmetry_array[0]=cplx_10(); - if ((2 == nop) || (4 == nop) || (8 == nop)) - return; - int ka = nop / 2; - -// int k_power = 2; -// for( ka = 1; k_power != nop; ka++ ) -// k_power *= 2; - + /* Plane symmetry admits two, four, or eight sections; the doubling below + cannot span the matrix for any other count. */ + if ((2 != nop) && (4 != nop) && (8 != nop)) + nec_stop("SYMMETRY ERROR - NOP: %d", nop ); + + /* Each pass doubles the width of the filled block, so the pass count is + log2(nop). */ + int ka = 1; + for (int k_power = 2; k_power != nop; k_power *= 2) + ka++; + for(int k = 0; k < ka; k++ ) { for(int i = 0; i < kk; i++ ) { for(int j = 0; j < kk; j++ ) { From 7475b206f4763592fbdcae96d0cf0b41860d3eb6 Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Sat, 1 Aug 2026 17:28:24 -0700 Subject: [PATCH 3/3] fix: reflect three symmetry planes in sequential Z, Y, X order reflect_plane() in c_geometry.cpp sized its copy machinery for two axes, using two axis-mask slots and a four-element tag table under a comment reading "2 for single axis, 4 for two axes", while a three-plane GX card sets num_copies to 8. The copy-axis selection chain covered only copies 1 through 3, so copies 4 through 7 kept their false-initialized flip flags and were written as verbatim duplicates of the original, and the tag table was indexed past its fourth element. A separate defect excluded the X axis from the tag increment shift, though every plane doubles the structure and so must double the increment. A three-plane deck therefore produced four distinct wires and four coincident copies of the first, with duplicate tags wired by the junction search into a chain the deck never described. - replace the two-axis mask/table pair with a copy_mask array built by iterating the requested axes in Z, Y, X order, doubling the filled block each pass, matching the sequential mirroring nec2c performs - derive flip_x/flip_y/flip_z for each copy from copy_mask bits instead of the four-case axis1_mask/axis2_mask chain, so every copy index up to 8 resolves correctly - derive tag_increment as itx * num_copies so X reflections double the tag increment along with Y and Z - replace the itagi + tag_inc[copy] lookup with itagi + copy * itx, removing the four-element table copy indices could exceed Signed-off-by: Eric Wheeler --- src/c_geometry.cpp | 75 ++++++++++++++-------------------------------- 1 file changed, 23 insertions(+), 52 deletions(-) diff --git a/src/c_geometry.cpp b/src/c_geometry.cpp index 1a68af1..f334ec6 100644 --- a/src/c_geometry.cpp +++ b/src/c_geometry.cpp @@ -1152,29 +1152,28 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment ) /* --- planar symmetry --- */ /* sym_plane bits: 0 = negate X, 1 = negate Y, 2 = negate Z */ - int mask = sym_plane; - bool neg_x = mask & 1; - bool neg_y = mask & 2; - bool neg_z = mask & 4; - int num_axes = (neg_x?1:0) + (neg_y?1:0) + (neg_z?1:0); - int num_copies = 1 << num_axes; /* 2 for single axis, 4 for two axes */ - int64_t orig_n = n_segments; int64_t orig_m = m; int itx = tag_increment; - /* Axis priority: Z (bit 2) > Y (bit 1) > X (bit 0). - For two-axis combos, copy order is: original, axis1-neg, axis2-neg, both-neg. - Tag increments: +0, +itx, +2*itx, +3*itx. */ - int axis1_mask = 0, axis2_mask = 0; - if ( neg_z ) { axis1_mask = 4; axis2_mask = neg_y ? 2 : (neg_x ? 1 : 0); } - else if ( neg_y ) { axis1_mask = 2; axis2_mask = neg_x ? 1 : 0; } - else if ( neg_x ) { axis1_mask = 1; } + /* Build copy masks in the Z, Y, X order used by sequential NEC-2 + reflection passes. Each pass appends a reflected copy of every existing + combination. */ + int copy_mask[8] = { 0 }; + int num_copies = 1; + + for ( int axis_mask = 4; axis_mask > 0; axis_mask >>= 1 ) + { + if ( 0 == (sym_plane & axis_mask) ) + continue; - int tag_inc[4] = { 0, itx, 2*itx, 3*itx }; + for ( int copy = 0; copy < num_copies; copy++ ) + copy_mask[num_copies + copy] = copy_mask[copy] | axis_mask; + + num_copies *= 2; + } - /* Final tag_increment: X never doubles, Z and Y always double */ - tag_increment = itx << ( (neg_z ? 1 : 0) + (neg_y ? 1 : 0) ); + tag_increment = itx * num_copies; /* --- SEGMENTS --- */ if ( orig_n > 0 ) @@ -1196,23 +1195,9 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment ) for( int copy = 1; copy < num_copies; copy++ ) { /* Determine which axes to negate for this copy */ - bool flip_x = false, flip_y = false, flip_z = false; - if ( copy == 1 ) - { - if ( axis1_mask == 4 ) flip_z = true; - else if ( axis1_mask == 2 ) flip_y = true; - else if ( axis1_mask == 1 ) flip_x = true; - } - else if ( copy == 2 ) - { - if ( axis2_mask == 4 ) flip_z = true; - else if ( axis2_mask == 2 ) flip_y = true; - else if ( axis2_mask == 1 ) flip_x = true; - } - else if ( copy == 3 ) - { - flip_x = neg_x; flip_y = neg_y; flip_z = neg_z; - } + bool flip_x = (0 != (copy_mask[copy] & 1)); + bool flip_y = (0 != (copy_mask[copy] & 2)); + bool flip_z = (0 != (copy_mask[copy] & 4)); int64_t base = copy * orig_n; @@ -1269,7 +1254,7 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment ) if ( itagi == 0 ) segment_tags[nx] = 0; if ( itagi != 0 ) - segment_tags[nx] = itagi + tag_inc[copy]; + segment_tags[nx] = itagi + copy * itx; segment_radius[nx] = segment_radius[i]; } @@ -1298,23 +1283,9 @@ void c_geometry::reflect_plane( int sym_plane, int& tag_increment ) for( int copy = 1; copy < num_copies; copy++ ) { - bool flip_x = false, flip_y = false, flip_z = false; - if ( copy == 1 ) - { - if ( axis1_mask == 4 ) flip_z = true; - else if ( axis1_mask == 2 ) flip_y = true; - else if ( axis1_mask == 1 ) flip_x = true; - } - else if ( copy == 2 ) - { - if ( axis2_mask == 4 ) flip_z = true; - else if ( axis2_mask == 2 ) flip_y = true; - else if ( axis2_mask == 1 ) flip_x = true; - } - else if ( copy == 3 ) - { - flip_x = neg_x; flip_y = neg_y; flip_z = neg_z; - } + bool flip_x = (0 != (copy_mask[copy] & 1)); + bool flip_y = (0 != (copy_mask[copy] & 2)); + bool flip_z = (0 != (copy_mask[copy] & 4)); int neg_count = (flip_x?1:0) + (flip_y?1:0) + (flip_z?1:0); bool neg_psalp = (neg_count % 2 == 1);