Skip to content

Divide by a subnormal pivot in banded LU factorization instead of scaling by its reciprocal (Reference-LAPACK PR 1379) #6177

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
martin-frbg merged 3 commits into from
Oct 26, 2026
+683 −341
Merged

Divide by a subnormal pivot in banded LU factorization instead of scaling by its reciprocal (Reference-LAPACK PR 1379) #6177

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
48c7b63
Divide by a subnormal pivot instead of scaling by its reciprocal (Ref…
martin-frbg Oct 26, 2026
b18e282
Divide by a subnormal pivot in banded LU factorization instead of sca…
martin-frbg Oct 26, 2026
6c35ade
Divide by a subnormal pivot in banded LU factorization instead of sca…
martin-frbg Oct 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 27 additions & 29 deletions lapack-netlib/SRC/cgbtf2.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,20 +249,6 @@ typedef struct Namelist Namelist;
*/


/* -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/



/* -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/



/* Table of constant values */

static complex c_b1 = {1.f,0.f};
Expand All @@ -276,7 +262,6 @@ e algorithm. */
/* Online html documentation available at */
/* http://www.netlib.org/lapack/explore-html/ */

/* > \htmlonly */
/* > Download CGBTF2 + dependencies */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbtf2.
f"> */
Expand All @@ -287,7 +272,6 @@ f"> */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbtf2.
f"> */
/* > [TXT]</a> */
/* > \endhtmlonly */

/* Definition: */
/* =========== */
Expand Down Expand Up @@ -385,9 +369,7 @@ f"> */
/* > \author Univ. of Colorado Denver */
/* > \author NAG Ltd. */

/* > \date December 2016 */

/* > \ingroup complexGBcomputational */
/* > \ingroup gbtf2 */

/* > \par Further Details: */
/* ===================== */
Expand Down Expand Up @@ -418,23 +400,24 @@ f"> */
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4;
complex q__1={0.,0.};
complex q__1;

/* Local variables */
integer i__, j;
integer i__, j, km, jp, ju, kv;
extern /* Subroutine */ void cscal_(integer *, complex *, complex *,
integer *), cgeru_(integer *, integer *, complex *, complex *,
integer *, complex *, integer *, complex *, integer *), cswap_(
integer *, complex *, integer *, complex *, integer *);
integer km, jp, ju, kv;
real sfmin;
extern /* Subroutine */ void cswap_(integer *, complex *, integer *,
complex *, integer *);
extern integer icamax_(integer *, complex *, integer *);
extern real slamch_(char *);
extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen);


/* -- LAPACK computational routine (version 3.7.0) -- */
/* -- LAPACK computational routine -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* December 2016 */


/* ===================================================================== */
Expand All @@ -445,7 +428,7 @@ f"> */

/* Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1 * 1;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
--ipiv;

Expand Down Expand Up @@ -478,6 +461,10 @@ f"> */
return;
}

/* Compute machine safe minimum */

sfmin = slamch_("S");

/* Gaussian elimination with partial pivoting */

/* Set fill-in elements in columns KU+2 to KV to zero. */
Expand Down Expand Up @@ -542,14 +529,25 @@ f"> */

/* Compute multipliers. */

c_div(&q__1, &c_b1, &ab[kv + 1 + j * ab_dim1]);
cscal_(&km, &q__1, &ab[kv + 2 + j * ab_dim1], &c__1);
if (c_abs(&ab[kv + 1 + j * ab_dim1]) >= sfmin) {
c_div(&q__1, &c_b1, &ab[kv + 1 + j * ab_dim1]);
cscal_(&km, &q__1, &ab[kv + 2 + j * ab_dim1], &c__1);
} else {
i__2 = km;
for (i__ = 1; i__ <= i__2; ++i__) {
i__3 = kv + 1 + i__ + j * ab_dim1;
c_div(&q__1, &ab[kv + 1 + i__ + j * ab_dim1], &ab[kv
+ 1 + j * ab_dim1]);
ab[i__3].r = q__1.r, ab[i__3].i = q__1.i;
/* L35: */
}
}

/* Update trailing submatrix within the band. */

if (ju > j) {
i__2 = ju - j;
q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__3 = *ldab - 1;
i__4 = *ldab - 1;
cgeru_(&km, &i__2, &q__1, &ab[kv + 2 + j * ab_dim1], &
Expand Down
22 changes: 17 additions & 5 deletions lapack-netlib/SRC/cgbtf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,13 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CGBTF2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbtf2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbtf2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbtf2.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -115,7 +113,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complexGBcomputational
*> \ingroup gbtf2
*
*> \par Further Details:
* =====================
Expand All @@ -142,6 +140,7 @@
*>
* =====================================================================
SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
IMPLICIT NONE
*
* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
Expand All @@ -163,11 +162,13 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
$ ZERO = ( 0.0E+0, 0.0E+0 ) )
* ..
* .. Local Scalars ..
REAL SFMIN
INTEGER I, J, JP, JU, KM, KV
* ..
* .. External Functions ..
REAL SLAMCH
INTEGER ICAMAX
EXTERNAL ICAMAX
EXTERNAL SLAMCH, ICAMAX
* ..
* .. External Subroutines ..
EXTERNAL CGERU, CSCAL, CSWAP, XERBLA
Expand Down Expand Up @@ -206,6 +207,10 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
IF( M.EQ.0 .OR. N.EQ.0 )
$ RETURN
*
* Compute machine safe minimum
*
SFMIN = SLAMCH('S')
*
* Gaussian elimination with partial pivoting
*
* Set fill-in elements in columns KU+2 to KV to zero.
Expand Down Expand Up @@ -249,7 +254,14 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
*
* Compute multipliers.
*
CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN
CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ),
$ 1 )
ELSE
DO 35 I = 1, KM
AB( KV+1+I, J ) = AB( KV+1+I, J ) / AB( KV+1, J )
35 CONTINUE
END IF
*
* Update trailing submatrix within the band.
*
Expand Down
71 changes: 34 additions & 37 deletions lapack-netlib/SRC/cgbtrf.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,20 +249,6 @@ typedef struct Namelist Namelist;
*/


/* -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/



/* -- translated by f2c (version 20000121).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/



/* Table of constant values */

static complex c_b1 = {1.f,0.f};
Expand All @@ -276,7 +262,6 @@ static integer c__65 = 65;
/* Online html documentation available at */
/* http://www.netlib.org/lapack/explore-html/ */

/* > \htmlonly */
/* > Download CGBTRF + dependencies */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbtrf.
f"> */
Expand All @@ -287,7 +272,6 @@ f"> */
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbtrf.
f"> */
/* > [TXT]</a> */
/* > \endhtmlonly */

/* Definition: */
/* =========== */
Expand Down Expand Up @@ -385,9 +369,7 @@ f"> */
/* > \author Univ. of Colorado Denver */
/* > \author NAG Ltd. */

/* > \date December 2016 */

/* > \ingroup complexGBcomputational */
/* > \ingroup gbtrf */

/* > \par Further Details: */
/* ===================== */
Expand Down Expand Up @@ -417,41 +399,41 @@ f"> */
{
/* System generated locals */
integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4, i__5, i__6;
complex q__1={0.,0.};
complex q__1;

/* Local variables */
integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju,
kv, nw;
complex temp;
integer i__, j;
extern /* Subroutine */ void cscal_(integer *, complex *, complex *,
integer *), cgemm_(char *, char *, integer *, integer *, integer *
, complex *, complex *, integer *, complex *, integer *, complex *
, complex *, integer *), cgeru_(integer *,
integer *, complex *, complex *, integer *, complex *, integer *,
complex *, integer *), ccopy_(integer *, complex *, integer *,
complex *, integer *);
real sfmin;
extern /* Subroutine */ void ccopy_(integer *, complex *, integer *,
complex *, integer *), cswap_(integer *, complex *, integer *,
complex *, integer *);
complex work13[4160] /* was [65][64] */, work31[4160] /*
was [65][64] */;
extern /* Subroutine */ void ctrsm_(char *, char *, char *, char *,
integer *, integer *, complex *, complex *, integer *, complex *,
integer *), cgbtf2_(integer *,
integer *, integer *, integer *, complex *, integer *, integer *,
integer *);
integer i2, i3, j2, j3, k2;
extern /* Subroutine */ void cgbtf2_(integer *, integer *, integer *,
integer *, complex *, integer *, integer *, integer *);
integer jb, nb, ii, jj, jm, ip, jp, km, ju, kv;
extern integer icamax_(integer *, complex *, integer *);
integer nw;
extern real slamch_(char *);
extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ void claswp_(integer *, complex *, integer *,
integer *, integer *, integer *, integer *);


/* -- LAPACK computational routine (version 3.7.0) -- */
/* -- LAPACK computational routine -- */
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
/* December 2016 */


/* ===================================================================== */
Expand All @@ -462,7 +444,7 @@ f"> */

/* Parameter adjustments */
ab_dim1 = *ldab;
ab_offset = 1 + ab_dim1 * 1;
ab_offset = 1 + ab_dim1;
ab -= ab_offset;
--ipiv;

Expand Down Expand Up @@ -495,6 +477,10 @@ f"> */
return;
}

/* Compute machine safe minimum */

sfmin = slamch_("S");

/* Determine the block size for this environment */

nb = ilaenv_(&c__1, "CGBTRF", " ", m, n, kl, ku, (ftnlen)6, (ftnlen)1);
Expand Down Expand Up @@ -650,8 +636,19 @@ f"> */

/* Compute multipliers */

c_div(&q__1, &c_b1, &ab[kv + 1 + jj * ab_dim1]);
cscal_(&km, &q__1, &ab[kv + 2 + jj * ab_dim1], &c__1);
if (c_abs(&ab[kv + 1 + jj * ab_dim1]) >= sfmin) {
c_div(&q__1, &c_b1, &ab[kv + 1 + jj * ab_dim1]);
cscal_(&km, &q__1, &ab[kv + 2 + jj * ab_dim1], &c__1);
} else {
i__4 = km;
for (i__ = 1; i__ <= i__4; ++i__) {
i__5 = kv + 1 + i__ + jj * ab_dim1;
c_div(&q__1, &ab[kv + 1 + i__ + jj * ab_dim1], &
ab[kv + 1 + jj * ab_dim1]);
ab[i__5].r = q__1.r, ab[i__5].i = q__1.i;
/* L75: */
}
}

/* Update trailing submatrix within the band and within */
/* the current block. JM is the index of the last column */
Expand All @@ -662,7 +659,7 @@ f"> */
jm = f2cmin(i__4,i__5);
if (jm > jj) {
i__4 = jm - jj;
q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__5 = *ldab - 1;
i__6 = *ldab - 1;
cgeru_(&km, &i__4, &q__1, &ab[kv + 2 + jj * ab_dim1],
Expand Down Expand Up @@ -756,7 +753,7 @@ f"> */

/* Update A22 */

q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__3 = *ldab - 1;
i__4 = *ldab - 1;
i__5 = *ldab - 1;
Expand All @@ -771,7 +768,7 @@ f"> */

/* Update A32 */

q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__3 = *ldab - 1;
i__4 = *ldab - 1;
cgemm_("No transpose", "No transpose", &i3, &j2, &jb,
Expand Down Expand Up @@ -810,7 +807,7 @@ f"> */

/* Update A23 */

q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__3 = *ldab - 1;
i__4 = *ldab - 1;
cgemm_("No transpose", "No transpose", &i2, &j3, &jb,
Expand All @@ -823,7 +820,7 @@ f"> */

/* Update A33 */

q__1.r = -1.f, q__1.i = 0.f;
q__1.r = -1.f, q__1.i = -0.f;
i__3 = *ldab - 1;
cgemm_("No transpose", "No transpose", &i3, &j3, &jb,
&q__1, work31, &c__65, work13, &c__65, &c_b1,
Expand Down
Loading
Loading