NF01BS DLACPY Overlapping Memory Bug

[jamestjsp]

Summary

The SLICOT routine NF01BS (QR factorization of Jacobian in compressed form) contains a bug where DLACPY is called with overlapping source and destination memory regions. This is undefined behavior per the LAPACK specification.

Bug Location

File: SLICOT-Reference/src/NF01BS.f
Line: 519

CALL DLACPY( 'Full', N, ST, J(LDJ*BSN+1), LDJ, J(IBSN), N )

Where IBSN = N*BSN + 1 (set at line 518).

Root Cause

The routine reshapes matrix J from leading dimension LDJ (= M = BNBSM) to leading dimension N (= BNBSN + ST). When LDJ > N, the source and destination regions overlap in memory.

Memory Layout

Source:      J(LDJ*BSN+1) with leading dimension LDJ
Destination: J(N*BSN+1)   with leading dimension N

When LDJ > N:
  Source offset  = LDJ * BSN = larger
  Dest offset    = N * BSN   = smaller

  These regions overlap!

Trigger Conditions

1. BN > 1 (multiple blocks - general case, not special case)
2. BSN > 0 (nonlinear part exists)
3. BSM > BSN (more rows than columns per block, MMN > 0)
4. ST > 0 (linear part exists)

LAPACK DLACPY Specification

From LAPACK documentation:

DLACPY copies all or part of a two-dimensional matrix A to another matrix B.

DLACPY does not handle overlapping memory regions. It uses simple nested loops (effectively memcpy semantics), not memmove.

Reproduction

Prerequisites

• gfortran
• LAPACK/BLAS libraries
• Valgrind

Build

cd tests/fortran
gfortran -g -O0 -fcheck=all -Wall -c \
    ../../SLICOT-Reference/src/NF01BS.f \
    ../../SLICOT-Reference/src/MD03BX.f
gfortran -g -O0 -fcheck=all -Wall -o test_nf01bs_overlap \
    test_nf01bs_overlap.f NF01BS.o MD03BX.o -llapack -lblas

Run with Valgrind

valgrind --track-origins=yes ./test_nf01bs_overlap

Expected Output

==19687== Source and destination overlap in memcpy(0x..., 0x..., 48)
==19687==    at 0x...: __GI_memcpy (...)
==19687==    by 0x...: dlacpy_ (liblapack.so)
==19687==    by 0x...: nf01bs_ (NF01BS.f:479)
==19687==    by 0x...: MAIN__ (test_nf01bs_overlap.f:96)

Test Parameters

The test uses parameters that reliably trigger the overlap:

Parameter	Value	Description
BN	2	Number of blocks
BSM	4	Block rows
BSN	2	Block columns
ST	2	Linear part size
N	6	Total columns (BN*BSN + ST)
M	8	Total rows (BN*BSM)
LDJ	8	Leading dimension (= M)

Memory offsets:

• Source: LDJ*BSN + 1 = 17 (1-based)
• Destination: N*BSN + 1 = 13 (1-based)
• Overlap because source starts at 17, destination at 13, with LDJ=8 > N=6

Fix in C Translation

The C translation (src/NF/nf01bs.c) uses a local helper function with memmove:

static void local_dlacpy_safe(i32 m, i32 n, const f64 *a, i32 lda, f64 *b, i32 ldb) {
    for (i32 j = 0; j < n; j++) {
        memmove(&b[j * ldb], &a[j * lda], m * sizeof(f64));
    }
}

This correctly handles overlapping memory regions.

Impact

• Severity: Medium
• Behavior: Undefined (may produce incorrect results silently)
• Detection: Only visible with memory debuggers (Valgrind, ASAN)
• Affected versions: All SLICOT versions containing NF01BS

Why It May Not Manifest

1. LAPACK's column-by-column copy order may accidentally work for some overlap patterns
2. Specific parameter combinations may avoid the overlap entirely
3. Resulting numerical errors may be small enough to go unnoticed

Recommendation

The original SLICOT Fortran code should be patched to use a custom copy routine that handles overlapping memory, similar to the C translation fix.

[jamestjsp]

C
C     Test program to reproduce DLACPY overlap bug in NF01BS
C
C     Bug location: NF01BS.f line 519
C     CALL DLACPY( 'Full', N, ST, J(LDJ*BSN+1), LDJ, J(IBSN), N )
C
C     The source and destination regions overlap when LDJ > N
C
      PROGRAM TEST_NF01BS_OVERLAP
      IMPLICIT NONE
C
C     Parameters to trigger overlap:
C     BN = 2, BSM = 4, BSN = 2, ST = 2
C     N = BN*BSN + ST = 6
C     M = BN*BSM = 8
C     LDJ = M = 8
C
      INTEGER BN, BSM, BSN, ST, N, M, LDJ_INIT
      PARAMETER ( BN = 2, BSM = 4, BSN = 2, ST = 2 )
      PARAMETER ( N = BN*BSN + ST )
      PARAMETER ( M = BN*BSM )
      PARAMETER ( LDJ_INIT = M )
C
C     Array dimensions
      INTEGER NC
      PARAMETER ( NC = BSN + ST )
C
C     Local variables
      INTEGER INFO, LDWORK, LIPAR, I, K, LDJ
      DOUBLE PRECISION FNORM, GNORM
C
C     Arrays
      INTEGER IPAR(4), IPVT(N)
      DOUBLE PRECISION J(LDJ_INIT, NC), J_ORIG(LDJ_INIT, NC)
      DOUBLE PRECISION E(M), JNORMS(N), DWORK(1000)
C
C     Initialize LDJ (must be variable since NF01BS modifies it)
      LDJ = LDJ_INIT
C
C     Initialize IPAR
      IPAR(1) = ST
      IPAR(2) = BN
      IPAR(3) = BSM
      IPAR(4) = BSN
      LIPAR = 4
C
C     Initialize J matrix with recognizable pattern
C     Column-major: J(i,j) stored at position i + (j-1)*LDJ
      DO K = 1, NC
         DO I = 1, M
            J(I,K) = DBLE(I + (K-1)*100)
            J_ORIG(I,K) = J(I,K)
         END DO
      END DO
C
C     Initialize error vector
      DO I = 1, M
         E(I) = 1.0D0
      END DO
      FNORM = SQRT(DBLE(M))
C
C     Workspace
      LDWORK = 1000
C
C     Print test setup
      WRITE(*,*) '=================================================='
      WRITE(*,*) 'NF01BS DLACPY Overlap Bug Reproduction Test'
      WRITE(*,*) '=================================================='
      WRITE(*,*)
      WRITE(*,*) 'Parameters:'
      WRITE(*,*) '  BN  =', BN,  ' (number of blocks)'
      WRITE(*,*) '  BSM =', BSM, ' (block rows)'
      WRITE(*,*) '  BSN =', BSN, ' (block cols)'
      WRITE(*,*) '  ST  =', ST,  ' (linear part size)'
      WRITE(*,*) '  N   =', N,   ' (total columns)'
      WRITE(*,*) '  M   =', M,   ' (total rows)'
      WRITE(*,*) '  LDJ =', LDJ, ' (leading dimension)'
      WRITE(*,*)
      WRITE(*,*) 'Memory layout analysis:'
      WRITE(*,*) '  Source offset  = LDJ*BSN+1 =', LDJ*BSN+1
      WRITE(*,*) '  Dest offset    = N*BSN+1   =', N*BSN+1
      WRITE(*,*) '  LDJ > N:', LDJ, '>', N, '=> OVERLAP!'
      WRITE(*,*)
C
C     Print original J matrix (last block column - linear part)
      WRITE(*,*) 'Original J matrix (cols BSN+1 to NC, linear part):'
      DO I = 1, M
         WRITE(*,'(A,I2,A,10F8.1)') '  Row ', I, ':',
     $        (J_ORIG(I,K), K=BSN+1, NC)
      END DO
      WRITE(*,*)
C
C     Call NF01BS
      WRITE(*,*) 'Calling NF01BS...'
      CALL NF01BS( N, IPAR, LIPAR, FNORM, J, LDJ, E, JNORMS,
     $             GNORM, IPVT, DWORK, LDWORK, INFO )
C
      WRITE(*,*)
      WRITE(*,*) 'Results:'
      WRITE(*,*) '  INFO  =', INFO
      WRITE(*,*) '  GNORM =', GNORM
      WRITE(*,*) '  LDJ (after) =', LDJ
      WRITE(*,*)
C
C     Print pivot indices
      WRITE(*,*) 'Pivot indices:'
      WRITE(*,'(A,20I4)') '  IPVT:', (IPVT(I), I=1,N)
      WRITE(*,*)
C
C     Print column norms
      WRITE(*,*) 'Column norms:'
      WRITE(*,'(A,6F10.4)') '  JNORMS:', (JNORMS(I), I=1,N)
      WRITE(*,*)
C
      WRITE(*,*) '=================================================='
      WRITE(*,*) 'Test completed. Run with valgrind to see overlap:'
      WRITE(*,*) '  valgrind --track-origins=yes ./test_nf01bs_overlap'
      WRITE(*,*) '=================================================='
C
      END

[jamestjsp]

SLICOT team may have a different opinion but in my experience the generative AI is a very useful tool especially to work with legacy code like this.

This and the previous bug I reported was identified by Claude Sonnet 4.5 using Claude Code.

[grisuthedragon]

I can confirm the bug. Obviously, this is a missuse of lacpy in this position and it only works, if it is a lacpy from reference lapack without any optimizations. Thus, I have to think of a clever replacement. Since aliasing in a fortran subroutine call is also not allowed (although it does not warn or throw an error)