From 4c6ffa55b3aece129e0ec3f4f611ec155d35e12b Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Sat, 13 Jun 2026 13:39:33 +0200 Subject: [PATCH] *trsm: remove zero-RHS skip for consistent NaN on singular matrix The left-solving non-transposed path in DTRSM/S/CTRSM/ZTRSM has an optimization that skips the division B(K,J)/A(K,K) when B(K,J) is exactly zero. When A(K,K) is also zero (singular matrix), this leaves B(K,J)=0 instead of computing 0/0=NaN, producing inconsistent results across different calling variants of the same triangular solve. Closes #636 --- BLAS/SRC/ctrsm.f | 36 ++++++++++++++++-------------------- BLAS/SRC/dtrsm.f | 24 ++++++++++-------------- BLAS/SRC/strsm.f | 36 ++++++++++++++++-------------------- BLAS/SRC/ztrsm.f | 36 ++++++++++++++++-------------------- 4 files changed, 58 insertions(+), 74 deletions(-) diff --git a/BLAS/SRC/ctrsm.f b/BLAS/SRC/ctrsm.f index 6275d4a59..39dececbb 100644 --- a/BLAS/SRC/ctrsm.f +++ b/BLAS/SRC/ctrsm.f @@ -283,29 +283,25 @@ SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B(I,J) = ALPHA*B(I,J) 30 CONTINUE END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF + DO 50 K = M,1,-1 + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 40 I = 1,K - 1 + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 40 CONTINUE 50 CONTINUE 60 CONTINUE - ELSE - DO 100 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 70 I = 1,M - B(I,J) = ALPHA*B(I,J) + ELSE + DO 100 J = 1,N + IF (ALPHA.NE.ONE) THEN + DO 70 I = 1,M + B(I,J) = ALPHA*B(I,J) 70 CONTINUE - END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF + END IF + DO 90 K = 1,M + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 80 I = K + 1,M + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 80 CONTINUE 90 CONTINUE 100 CONTINUE END IF diff --git a/BLAS/SRC/dtrsm.f b/BLAS/SRC/dtrsm.f index 0cd4a91a2..4fb64bf66 100644 --- a/BLAS/SRC/dtrsm.f +++ b/BLAS/SRC/dtrsm.f @@ -281,13 +281,11 @@ SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B(I,J) = ALPHA*B(I,J) 30 CONTINUE END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF + DO 50 K = M,1,-1 + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 40 I = 1,K - 1 + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 40 CONTINUE 50 CONTINUE 60 CONTINUE ELSE @@ -297,13 +295,11 @@ SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B(I,J) = ALPHA*B(I,J) 70 CONTINUE END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF + DO 90 K = 1,M + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 80 I = K + 1,M + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 80 CONTINUE 90 CONTINUE 100 CONTINUE END IF diff --git a/BLAS/SRC/strsm.f b/BLAS/SRC/strsm.f index 45fbb315b..495ee6fa3 100644 --- a/BLAS/SRC/strsm.f +++ b/BLAS/SRC/strsm.f @@ -281,29 +281,25 @@ SUBROUTINE STRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B(I,J) = ALPHA*B(I,J) 30 CONTINUE END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF + DO 50 K = M,1,-1 + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 40 I = 1,K - 1 + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 40 CONTINUE 50 CONTINUE 60 CONTINUE - ELSE - DO 100 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 70 I = 1,M - B(I,J) = ALPHA*B(I,J) + ELSE + DO 100 J = 1,N + IF (ALPHA.NE.ONE) THEN + DO 70 I = 1,M + B(I,J) = ALPHA*B(I,J) 70 CONTINUE - END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF + END IF + DO 90 K = 1,M + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 80 I = K + 1,M + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 80 CONTINUE 90 CONTINUE 100 CONTINUE END IF diff --git a/BLAS/SRC/ztrsm.f b/BLAS/SRC/ztrsm.f index 192190f6c..241fbcb30 100644 --- a/BLAS/SRC/ztrsm.f +++ b/BLAS/SRC/ztrsm.f @@ -283,29 +283,25 @@ SUBROUTINE ZTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B(I,J) = ALPHA*B(I,J) 30 CONTINUE END IF - DO 50 K = M,1,-1 - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 40 I = 1,K - 1 - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 40 CONTINUE - END IF + DO 50 K = M,1,-1 + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 40 I = 1,K - 1 + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 40 CONTINUE 50 CONTINUE 60 CONTINUE - ELSE - DO 100 J = 1,N - IF (ALPHA.NE.ONE) THEN - DO 70 I = 1,M - B(I,J) = ALPHA*B(I,J) + ELSE + DO 100 J = 1,N + IF (ALPHA.NE.ONE) THEN + DO 70 I = 1,M + B(I,J) = ALPHA*B(I,J) 70 CONTINUE - END IF - DO 90 K = 1,M - IF (B(K,J).NE.ZERO) THEN - IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) - DO 80 I = K + 1,M - B(I,J) = B(I,J) - B(K,J)*A(I,K) - 80 CONTINUE - END IF + END IF + DO 90 K = 1,M + IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) + DO 80 I = K + 1,M + B(I,J) = B(I,J) - B(K,J)*A(I,K) + 80 CONTINUE 90 CONTINUE 100 CONTINUE END IF