diff --git a/lapack-netlib/SRC/cgbtf2.c b/lapack-netlib/SRC/cgbtf2.c index 6f791d034c..7b98ac157e 100644 --- a/lapack-netlib/SRC/cgbtf2.c +++ b/lapack-netlib/SRC/cgbtf2.c @@ -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}; @@ -276,7 +262,6 @@ e algorithm. */ /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download CGBTF2 + dependencies */ /* > */ @@ -287,7 +272,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -385,9 +369,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complexGBcomputational */ +/* > \ingroup gbtf2 */ /* > \par Further Details: */ /* ===================== */ @@ -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 */ /* ===================================================================== */ @@ -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; @@ -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. */ @@ -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], & diff --git a/lapack-netlib/SRC/cgbtf2.f b/lapack-netlib/SRC/cgbtf2.f index 8c58d6e118..fdcfdb178a 100644 --- a/lapack-netlib/SRC/cgbtf2.f +++ b/lapack-netlib/SRC/cgbtf2.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CGBTF2 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexGBcomputational +*> \ingroup gbtf2 * *> \par Further Details: * ===================== @@ -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, -- @@ -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 @@ -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. @@ -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. * diff --git a/lapack-netlib/SRC/cgbtrf.c b/lapack-netlib/SRC/cgbtrf.c index 64d74b0c3b..d2473e7576 100644 --- a/lapack-netlib/SRC/cgbtrf.c +++ b/lapack-netlib/SRC/cgbtrf.c @@ -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}; @@ -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 */ /* > */ @@ -287,7 +272,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -385,9 +369,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complexGBcomputational */ +/* > \ingroup gbtrf */ /* > \par Further Details: */ /* ===================== */ @@ -417,30 +399,31 @@ 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); @@ -448,10 +431,9 @@ f"> */ 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 */ /* ===================================================================== */ @@ -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; @@ -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); @@ -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 */ @@ -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], @@ -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; @@ -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, @@ -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, @@ -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, diff --git a/lapack-netlib/SRC/cgbtrf.f b/lapack-netlib/SRC/cgbtrf.f index 8533e0061d..eb00b7fb14 100644 --- a/lapack-netlib/SRC/cgbtrf.f +++ b/lapack-netlib/SRC/cgbtrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CGBTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexGBcomputational +*> \ingroup gbtrf * *> \par Further Details: * ===================== @@ -141,6 +139,7 @@ *> * ===================================================================== SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -167,17 +166,20 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW COMPLEX TEMP + REAL SFMIN * .. * .. Local Arrays .. COMPLEX WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + REAL SLAMCH INTEGER ICAMAX, ILAENV - EXTERNAL ICAMAX, ILAENV + EXTERNAL SLAMCH, ICAMAX, ILAENV * .. * .. External Subroutines .. - EXTERNAL CCOPY, CGBTF2, CGEMM, CGERU, CLASWP, CSCAL, + EXTERNAL CCOPY, CGBTF2, CGEMM, CGERU, CLASWP, + $ CSCAL, $ CSWAP, CTRSM, XERBLA * .. * .. Intrinsic Functions .. @@ -214,6 +216,10 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'CGBTRF', ' ', M, N, KL, KU ) @@ -325,8 +331,15 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL CSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL CSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column @@ -394,7 +407,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A12 * - CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL CTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1 ) * @@ -402,7 +416,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A22 * - CALL CGEMM( 'No transpose', 'No transpose', I2, J2, + CALL CGEMM( 'No transpose', 'No transpose', I2, + $ J2, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+1, J+JB ), LDAB-1 ) @@ -412,7 +427,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A32 * - CALL CGEMM( 'No transpose', 'No transpose', I3, J2, + CALL CGEMM( 'No transpose', 'No transpose', I3, + $ J2, $ JB, -ONE, WORK31, LDWORK, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) @@ -432,7 +448,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A13 in the work array * - CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL CTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, $ WORK13, LDWORK ) * @@ -440,7 +457,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A23 * - CALL CGEMM( 'No transpose', 'No transpose', I2, J3, + CALL CGEMM( 'No transpose', 'No transpose', I2, + $ J3, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), $ LDAB-1 ) @@ -450,7 +468,8 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A33 * - CALL CGEMM( 'No transpose', 'No transpose', I3, J3, + CALL CGEMM( 'No transpose', 'No transpose', I3, + $ J3, $ JB, -ONE, WORK31, LDWORK, WORK13, $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF diff --git a/lapack-netlib/SRC/dgbtf2.c b/lapack-netlib/SRC/dgbtf2.c index 0ad61a4976..cfd90c1d60 100644 --- a/lapack-netlib/SRC/dgbtf2.c +++ b/lapack-netlib/SRC/dgbtf2.c @@ -226,12 +226,10 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; -static doublereal c_b9 = -1.; +static doublereal c_b11 = -1.; /* > \brief \b DGBTF2 computes the LU factorization of a general band matrix using the unblocked version of th e algorithm. */ @@ -241,7 +239,6 @@ e algorithm. */ /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download DGBTF2 + dependencies */ /* > */ @@ -252,7 +249,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -350,9 +346,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup doubleGBcomputational */ +/* > \ingroup gbtf2 */ /* > \par Further Details: */ /* ===================== */ @@ -386,22 +380,22 @@ f"> */ doublereal d__1; /* Local variables */ + integer i__, j, km, jp, ju, kv; extern /* Subroutine */ void dger_(integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, - integer *); - integer i__, j; - extern /* Subroutine */ void dscal_(integer *, doublereal *, doublereal *, - integer *), dswap_(integer *, doublereal *, integer *, doublereal - *, integer *); - integer km, jp, ju, kv; + integer *), dscal_(integer *, doublereal *, doublereal *, integer + *); + doublereal sfmin; + extern /* Subroutine */ void dswap_(integer *, doublereal *, integer *, + doublereal *, integer *); + extern doublereal dlamch_(char *); extern integer idamax_(integer *, doublereal *, integer *); 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 */ /* ===================================================================== */ @@ -412,7 +406,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -445,6 +439,10 @@ f"> */ return; } +/* Compute machine safe minimum */ + + sfmin = dlamch_("S"); + /* Gaussian elimination with partial pivoting */ /* Set fill-in elements in columns KU+2 to KV to zero. */ @@ -507,8 +505,17 @@ f"> */ /* Compute multipliers. */ - d__1 = 1. / ab[kv + 1 + j * ab_dim1]; - dscal_(&km, &d__1, &ab[kv + 2 + j * ab_dim1], &c__1); + if ((d__1 = ab[kv + 1 + j * ab_dim1], abs(d__1)) >= sfmin) { + d__1 = 1. / ab[kv + 1 + j * ab_dim1]; + dscal_(&km, &d__1, &ab[kv + 2 + j * ab_dim1], &c__1); + } else { + i__2 = km; + for (i__ = 1; i__ <= i__2; ++i__) { + ab[kv + 1 + i__ + j * ab_dim1] /= ab[kv + 1 + j * + ab_dim1]; +/* L35: */ + } + } /* Update trailing submatrix within the band. */ @@ -516,9 +523,9 @@ f"> */ i__2 = ju - j; i__3 = *ldab - 1; i__4 = *ldab - 1; - dger_(&km, &i__2, &c_b9, &ab[kv + 2 + j * ab_dim1], &c__1, - &ab[kv + (j + 1) * ab_dim1], &i__3, &ab[kv + 1 + - (j + 1) * ab_dim1], &i__4); + dger_(&km, &i__2, &c_b11, &ab[kv + 2 + j * ab_dim1], & + c__1, &ab[kv + (j + 1) * ab_dim1], &i__3, &ab[kv + + 1 + (j + 1) * ab_dim1], &i__4); } } } else { diff --git a/lapack-netlib/SRC/dgbtf2.f b/lapack-netlib/SRC/dgbtf2.f index a643682582..6b1991e763 100644 --- a/lapack-netlib/SRC/dgbtf2.f +++ b/lapack-netlib/SRC/dgbtf2.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DGBTF2 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleGBcomputational +*> \ingroup gbtf2 * *> \par Further Details: * ===================== @@ -142,6 +140,7 @@ *> * ===================================================================== SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -162,11 +161,13 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. + DOUBLE PRECISION SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IDAMAX - EXTERNAL IDAMAX + EXTERNAL DLAMCH, IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA @@ -205,6 +206,10 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -249,7 +254,14 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL DSCAL( 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. * diff --git a/lapack-netlib/SRC/dgbtrf.c b/lapack-netlib/SRC/dgbtrf.c index 6b78cd746a..430e940f79 100644 --- a/lapack-netlib/SRC/dgbtrf.c +++ b/lapack-netlib/SRC/dgbtrf.c @@ -226,14 +226,12 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c__65 = 65; -static doublereal c_b18 = -1.; -static doublereal c_b31 = 1.; +static doublereal c_b20 = -1.; +static doublereal c_b33 = 1.; /* > \brief \b DGBTRF */ @@ -242,7 +240,6 @@ static doublereal c_b31 = 1.; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download DGBTRF + dependencies */ /* > */ @@ -253,7 +250,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -351,9 +347,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup doubleGBcomputational */ +/* > \ingroup gbtrf */ /* > \par Further Details: */ /* ===================== */ @@ -386,29 +380,29 @@ f"> */ doublereal d__1; /* Local variables */ + integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju, + kv, nw; extern /* Subroutine */ void dger_(integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *); doublereal temp; - integer i__, j; extern /* Subroutine */ void dscal_(integer *, doublereal *, doublereal *, integer *), dgemm_(char *, char *, integer *, integer *, integer * , doublereal *, doublereal *, integer *, doublereal *, integer *, - doublereal *, doublereal *, integer *), dcopy_( - integer *, doublereal *, integer *, doublereal *, integer *), - dswap_(integer *, doublereal *, integer *, doublereal *, integer * - ); + doublereal *, doublereal *, integer *); + doublereal sfmin; + extern /* Subroutine */ void dcopy_(integer *, doublereal *, integer *, + doublereal *, integer *), dswap_(integer *, doublereal *, integer + *, doublereal *, integer *); doublereal work13[4160] /* was [65][64] */, work31[4160] /* was [65][64] */; extern /* Subroutine */ void dtrsm_(char *, char *, char *, char *, integer *, integer *, doublereal *, doublereal *, integer *, - doublereal *, integer *); - integer i2, i3, j2, j3, k2; - extern /* Subroutine */ void dgbtf2_(integer *, integer *, integer *, - integer *, doublereal *, integer *, integer *, integer *); - integer jb, nb, ii, jj, jm, ip, jp, km, ju, kv; + doublereal *, integer *), dgbtf2_( + integer *, integer *, integer *, integer *, doublereal *, integer + *, integer *, integer *); + extern doublereal dlamch_(char *); extern integer idamax_(integer *, doublereal *, integer *); - integer nw; extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); @@ -416,10 +410,9 @@ f"> */ 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 */ /* ===================================================================== */ @@ -430,7 +423,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -463,6 +456,10 @@ f"> */ return; } +/* Compute machine safe minimum */ + + sfmin = dlamch_("S"); + /* Determine the block size for this environment */ nb = ilaenv_(&c__1, "DGBTRF", " ", m, n, kl, ku, (ftnlen)6, (ftnlen)1); @@ -613,8 +610,18 @@ f"> */ /* Compute multipliers */ - d__1 = 1. / ab[kv + 1 + jj * ab_dim1]; - dscal_(&km, &d__1, &ab[kv + 2 + jj * ab_dim1], &c__1); + if ((d__1 = ab[kv + 1 + jj * ab_dim1], abs(d__1)) >= + sfmin) { + d__1 = 1. / ab[kv + 1 + jj * ab_dim1]; + dscal_(&km, &d__1, &ab[kv + 2 + jj * ab_dim1], &c__1); + } else { + i__4 = km; + for (i__ = 1; i__ <= i__4; ++i__) { + ab[kv + 1 + i__ + jj * ab_dim1] /= ab[kv + 1 + jj + * ab_dim1]; +/* L75: */ + } + } /* Update trailing submatrix within the band and within */ /* the current block. JM is the index of the last column */ @@ -627,7 +634,7 @@ f"> */ i__4 = jm - jj; i__5 = *ldab - 1; i__6 = *ldab - 1; - dger_(&km, &i__4, &c_b18, &ab[kv + 2 + jj * ab_dim1], + dger_(&km, &i__4, &c_b20, &ab[kv + 2 + jj * ab_dim1], &c__1, &ab[kv + (jj + 1) * ab_dim1], &i__5, & ab[kv + 1 + (jj + 1) * ab_dim1], &i__6); } @@ -708,7 +715,7 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; dtrsm_("Left", "Lower", "No transpose", "Unit", &jb, &j2, - &c_b31, &ab[kv + 1 + j * ab_dim1], &i__3, &ab[kv + &c_b33, &ab[kv + 1 + j * ab_dim1], &i__3, &ab[kv + 1 - jb + (j + jb) * ab_dim1], &i__4); if (i2 > 0) { @@ -719,9 +726,9 @@ f"> */ i__4 = *ldab - 1; i__5 = *ldab - 1; dgemm_("No transpose", "No transpose", &i2, &j2, &jb, - &c_b18, &ab[kv + 1 + jb + j * ab_dim1], &i__3, + &c_b20, &ab[kv + 1 + jb + j * ab_dim1], &i__3, &ab[kv + 1 - jb + (j + jb) * ab_dim1], &i__4, - &c_b31, &ab[kv + 1 + (j + jb) * ab_dim1], & + &c_b33, &ab[kv + 1 + (j + jb) * ab_dim1], & i__5); } @@ -732,8 +739,8 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; dgemm_("No transpose", "No transpose", &i3, &j2, &jb, - &c_b18, work31, &c__65, &ab[kv + 1 - jb + (j - + jb) * ab_dim1], &i__3, &c_b31, &ab[kv + *kl + &c_b20, work31, &c__65, &ab[kv + 1 - jb + (j + + jb) * ab_dim1], &i__3, &c_b33, &ab[kv + *kl + 1 - jb + (j + jb) * ab_dim1], &i__4); } } @@ -758,7 +765,7 @@ f"> */ i__3 = *ldab - 1; dtrsm_("Left", "Lower", "No transpose", "Unit", &jb, &j3, - &c_b31, &ab[kv + 1 + j * ab_dim1], &i__3, work13, + &c_b33, &ab[kv + 1 + j * ab_dim1], &i__3, work13, &c__65); if (i2 > 0) { @@ -768,8 +775,8 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; dgemm_("No transpose", "No transpose", &i2, &j3, &jb, - &c_b18, &ab[kv + 1 + jb + j * ab_dim1], &i__3, - work13, &c__65, &c_b31, &ab[jb + 1 + (j + kv) + &c_b20, &ab[kv + 1 + jb + j * ab_dim1], &i__3, + work13, &c__65, &c_b33, &ab[jb + 1 + (j + kv) * ab_dim1], &i__4); } @@ -779,8 +786,8 @@ f"> */ i__3 = *ldab - 1; dgemm_("No transpose", "No transpose", &i3, &j3, &jb, - &c_b18, work31, &c__65, work13, &c__65, & - c_b31, &ab[*kl + 1 + (j + kv) * ab_dim1], & + &c_b20, work31, &c__65, work13, &c__65, & + c_b33, &ab[*kl + 1 + (j + kv) * ab_dim1], & i__3); } diff --git a/lapack-netlib/SRC/dgbtrf.f b/lapack-netlib/SRC/dgbtrf.f index af305e26b4..8bee3120ca 100644 --- a/lapack-netlib/SRC/dgbtrf.f +++ b/lapack-netlib/SRC/dgbtrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DGBTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleGBcomputational +*> \ingroup gbtrf * *> \par Further Details: * ===================== @@ -141,6 +139,7 @@ *> * ===================================================================== SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -165,18 +164,20 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW - DOUBLE PRECISION TEMP + DOUBLE PRECISION SFMIN, TEMP * .. * .. Local Arrays .. DOUBLE PRECISION WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IDAMAX, ILAENV - EXTERNAL IDAMAX, ILAENV + EXTERNAL DLAMCH, IDAMAX, ILAENV * .. * .. External Subroutines .. - EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, DSCAL, + EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, + $ DSCAL, $ DSWAP, DTRSM, XERBLA * .. * .. Intrinsic Functions .. @@ -213,6 +214,10 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU ) @@ -324,8 +329,15 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL DSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column @@ -393,7 +405,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A12 * - CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL DTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1 ) * @@ -401,7 +414,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A22 * - CALL DGEMM( 'No transpose', 'No transpose', I2, J2, + CALL DGEMM( 'No transpose', 'No transpose', I2, + $ J2, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+1, J+JB ), LDAB-1 ) @@ -411,7 +425,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A32 * - CALL DGEMM( 'No transpose', 'No transpose', I3, J2, + CALL DGEMM( 'No transpose', 'No transpose', I3, + $ J2, $ JB, -ONE, WORK31, LDWORK, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) @@ -431,7 +446,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A13 in the work array * - CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL DTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, $ WORK13, LDWORK ) * @@ -439,7 +455,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A23 * - CALL DGEMM( 'No transpose', 'No transpose', I2, J3, + CALL DGEMM( 'No transpose', 'No transpose', I2, + $ J3, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), $ LDAB-1 ) @@ -449,7 +466,8 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A33 * - CALL DGEMM( 'No transpose', 'No transpose', I3, J3, + CALL DGEMM( 'No transpose', 'No transpose', I3, + $ J3, $ JB, -ONE, WORK31, LDWORK, WORK13, $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF diff --git a/lapack-netlib/SRC/sgbtf2.c b/lapack-netlib/SRC/sgbtf2.c index eb89cbfb07..ced33a932a 100644 --- a/lapack-netlib/SRC/sgbtf2.c +++ b/lapack-netlib/SRC/sgbtf2.c @@ -219,12 +219,10 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; -static real c_b9 = -1.f; +static real c_b11 = -1.f; /* > \brief \b SGBTF2 computes the LU factorization of a general band matrix using the unblocked version of th e algorithm. */ @@ -234,7 +232,6 @@ e algorithm. */ /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download SGBTF2 + dependencies */ /* > */ @@ -245,7 +242,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -343,9 +339,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup realGBcomputational */ +/* > \ingroup gbtf2 */ /* > \par Further Details: */ /* ===================== */ @@ -379,20 +373,21 @@ f"> */ real r__1; /* Local variables */ + integer i__, j, km, jp, ju, kv; extern /* Subroutine */ void sger_(integer *, integer *, real *, real *, - integer *, real *, integer *, real *, integer *); - integer i__, j; - extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *), - sswap_(integer *, real *, integer *, real *, integer *); - integer km, jp, ju, kv; + integer *, real *, integer *, real *, integer *), sscal_(integer * + , real *, real *, integer *); + real sfmin; + extern /* Subroutine */ void sswap_(integer *, real *, integer *, real *, + integer *); + extern real slamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer isamax_(integer *, real *, 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 */ /* ===================================================================== */ @@ -403,7 +398,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -436,6 +431,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. */ @@ -498,8 +497,17 @@ f"> */ /* Compute multipliers. */ - r__1 = 1.f / ab[kv + 1 + j * ab_dim1]; - sscal_(&km, &r__1, &ab[kv + 2 + j * ab_dim1], &c__1); + if ((r__1 = ab[kv + 1 + j * ab_dim1], abs(r__1)) >= sfmin) { + r__1 = 1.f / ab[kv + 1 + j * ab_dim1]; + sscal_(&km, &r__1, &ab[kv + 2 + j * ab_dim1], &c__1); + } else { + i__2 = km; + for (i__ = 1; i__ <= i__2; ++i__) { + ab[kv + 1 + i__ + j * ab_dim1] /= ab[kv + 1 + j * + ab_dim1]; +/* L35: */ + } + } /* Update trailing submatrix within the band. */ @@ -507,9 +515,9 @@ f"> */ i__2 = ju - j; i__3 = *ldab - 1; i__4 = *ldab - 1; - sger_(&km, &i__2, &c_b9, &ab[kv + 2 + j * ab_dim1], &c__1, - &ab[kv + (j + 1) * ab_dim1], &i__3, &ab[kv + 1 + - (j + 1) * ab_dim1], &i__4); + sger_(&km, &i__2, &c_b11, &ab[kv + 2 + j * ab_dim1], & + c__1, &ab[kv + (j + 1) * ab_dim1], &i__3, &ab[kv + + 1 + (j + 1) * ab_dim1], &i__4); } } } else { diff --git a/lapack-netlib/SRC/sgbtf2.f b/lapack-netlib/SRC/sgbtf2.f index 8bdea3902f..53026dd266 100644 --- a/lapack-netlib/SRC/sgbtf2.f +++ b/lapack-netlib/SRC/sgbtf2.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SGBTF2 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realGBcomputational +*> \ingroup gbtf2 * *> \par Further Details: * ===================== @@ -142,6 +140,7 @@ *> * ===================================================================== SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -162,11 +161,13 @@ SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. + REAL SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + REAL SLAMCH INTEGER ISAMAX - EXTERNAL ISAMAX + EXTERNAL SLAMCH, ISAMAX * .. * .. External Subroutines .. EXTERNAL SGER, SSCAL, SSWAP, XERBLA @@ -205,6 +206,10 @@ SUBROUTINE SGBTF2( 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. @@ -249,7 +254,14 @@ SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL SSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL SSCAL( 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. * diff --git a/lapack-netlib/SRC/sgbtrf.c b/lapack-netlib/SRC/sgbtrf.c index 3b5c34b898..a8075b477f 100644 --- a/lapack-netlib/SRC/sgbtrf.c +++ b/lapack-netlib/SRC/sgbtrf.c @@ -39,9 +39,6 @@ typedef float real; typedef double doublereal; typedef struct { real r, i; } complex; typedef struct { doublereal r, i; } doublecomplex; -typedef blasint logical; - -typedef char logical1; typedef char integer1; #define TRUE_ (1) @@ -222,14 +219,12 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c__65 = 65; -static real c_b18 = -1.f; -static real c_b31 = 1.f; +static real c_b20 = -1.f; +static real c_b33 = 1.f; /* > \brief \b SGBTRF */ @@ -238,7 +233,6 @@ static real c_b31 = 1.f; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download SGBTRF + dependencies */ /* > */ @@ -249,7 +243,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -347,9 +340,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup realGBcomputational */ +/* > \ingroup gbtrf */ /* > \par Further Details: */ /* ===================== */ @@ -382,22 +373,22 @@ f"> */ real r__1; /* Local variables */ + integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju, + kv, nw; extern /* Subroutine */ void sger_(integer *, integer *, real *, real *, integer *, real *, integer *, real *, integer *); real temp; - integer i__, j; extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *), sgemm_(char *, char *, integer *, integer *, integer *, real *, real *, integer *, real *, integer *, real *, real *, integer *); - real work13[4160] /* was [65][64] */, work31[4160] /* was [65][ - 64] */; - integer i2, i3, j2, j3, k2; + real sfmin, work13[4160] /* was [65][64] */, work31[4160] /* + was [65][64] */; extern /* Subroutine */ void scopy_(integer *, real *, integer *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer * ), strsm_(char *, char *, char *, char *, integer *, integer *, real *, real *, integer *, real *, integer *), sgbtf2_(integer *, integer *, integer *, integer *, real *, integer *, integer *, integer *); - integer jb, nb, ii, jj, jm, ip, jp, km, ju, kv, nw; + extern real slamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen), isamax_(integer *, real *, @@ -406,10 +397,9 @@ f"> */ *, 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 */ /* ===================================================================== */ @@ -420,7 +410,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -453,6 +443,10 @@ f"> */ return; } +/* Compute machine safe minimum */ + + sfmin = slamch_("S"); + /* Determine the block size for this environment */ nb = ilaenv_(&c__1, "SGBTRF", " ", m, n, kl, ku, (ftnlen)6, (ftnlen)1); @@ -603,8 +597,18 @@ f"> */ /* Compute multipliers */ - r__1 = 1.f / ab[kv + 1 + jj * ab_dim1]; - sscal_(&km, &r__1, &ab[kv + 2 + jj * ab_dim1], &c__1); + if ((r__1 = ab[kv + 1 + jj * ab_dim1], abs(r__1)) >= + sfmin) { + r__1 = 1.f / ab[kv + 1 + jj * ab_dim1]; + sscal_(&km, &r__1, &ab[kv + 2 + jj * ab_dim1], &c__1); + } else { + i__4 = km; + for (i__ = 1; i__ <= i__4; ++i__) { + ab[kv + 1 + i__ + jj * ab_dim1] /= ab[kv + 1 + jj + * ab_dim1]; +/* L75: */ + } + } /* Update trailing submatrix within the band and within */ /* the current block. JM is the index of the last column */ @@ -617,7 +621,7 @@ f"> */ i__4 = jm - jj; i__5 = *ldab - 1; i__6 = *ldab - 1; - sger_(&km, &i__4, &c_b18, &ab[kv + 2 + jj * ab_dim1], + sger_(&km, &i__4, &c_b20, &ab[kv + 2 + jj * ab_dim1], &c__1, &ab[kv + (jj + 1) * ab_dim1], &i__5, & ab[kv + 1 + (jj + 1) * ab_dim1], &i__6); } @@ -698,7 +702,7 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; strsm_("Left", "Lower", "No transpose", "Unit", &jb, &j2, - &c_b31, &ab[kv + 1 + j * ab_dim1], &i__3, &ab[kv + &c_b33, &ab[kv + 1 + j * ab_dim1], &i__3, &ab[kv + 1 - jb + (j + jb) * ab_dim1], &i__4); if (i2 > 0) { @@ -709,9 +713,9 @@ f"> */ i__4 = *ldab - 1; i__5 = *ldab - 1; sgemm_("No transpose", "No transpose", &i2, &j2, &jb, - &c_b18, &ab[kv + 1 + jb + j * ab_dim1], &i__3, + &c_b20, &ab[kv + 1 + jb + j * ab_dim1], &i__3, &ab[kv + 1 - jb + (j + jb) * ab_dim1], &i__4, - &c_b31, &ab[kv + 1 + (j + jb) * ab_dim1], & + &c_b33, &ab[kv + 1 + (j + jb) * ab_dim1], & i__5); } @@ -722,8 +726,8 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; sgemm_("No transpose", "No transpose", &i3, &j2, &jb, - &c_b18, work31, &c__65, &ab[kv + 1 - jb + (j - + jb) * ab_dim1], &i__3, &c_b31, &ab[kv + *kl + &c_b20, work31, &c__65, &ab[kv + 1 - jb + (j + + jb) * ab_dim1], &i__3, &c_b33, &ab[kv + *kl + 1 - jb + (j + jb) * ab_dim1], &i__4); } } @@ -748,7 +752,7 @@ f"> */ i__3 = *ldab - 1; strsm_("Left", "Lower", "No transpose", "Unit", &jb, &j3, - &c_b31, &ab[kv + 1 + j * ab_dim1], &i__3, work13, + &c_b33, &ab[kv + 1 + j * ab_dim1], &i__3, work13, &c__65); if (i2 > 0) { @@ -758,8 +762,8 @@ f"> */ i__3 = *ldab - 1; i__4 = *ldab - 1; sgemm_("No transpose", "No transpose", &i2, &j3, &jb, - &c_b18, &ab[kv + 1 + jb + j * ab_dim1], &i__3, - work13, &c__65, &c_b31, &ab[jb + 1 + (j + kv) + &c_b20, &ab[kv + 1 + jb + j * ab_dim1], &i__3, + work13, &c__65, &c_b33, &ab[jb + 1 + (j + kv) * ab_dim1], &i__4); } @@ -769,8 +773,8 @@ f"> */ i__3 = *ldab - 1; sgemm_("No transpose", "No transpose", &i3, &j3, &jb, - &c_b18, work31, &c__65, work13, &c__65, & - c_b31, &ab[*kl + 1 + (j + kv) * ab_dim1], & + &c_b20, work31, &c__65, work13, &c__65, & + c_b33, &ab[*kl + 1 + (j + kv) * ab_dim1], & i__3); } diff --git a/lapack-netlib/SRC/sgbtrf.f b/lapack-netlib/SRC/sgbtrf.f index 56712f157e..c9ac5327b0 100644 --- a/lapack-netlib/SRC/sgbtrf.f +++ b/lapack-netlib/SRC/sgbtrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SGBTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realGBcomputational +*> \ingroup gbtrf * *> \par Further Details: * ===================== @@ -141,6 +139,7 @@ *> * ===================================================================== SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -165,18 +164,20 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW - REAL TEMP + REAL SFMIN, TEMP * .. * .. Local Arrays .. REAL WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + REAL SLAMCH INTEGER ILAENV, ISAMAX - EXTERNAL ILAENV, ISAMAX + EXTERNAL SLAMCH, ILAENV, ISAMAX * .. * .. External Subroutines .. - EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, SSCAL, + EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, + $ SSCAL, $ SSWAP, STRSM, XERBLA * .. * .. Intrinsic Functions .. @@ -213,6 +214,10 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU ) @@ -324,8 +329,15 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL SSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column @@ -393,7 +405,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A12 * - CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL STRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1 ) * @@ -401,7 +414,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A22 * - CALL SGEMM( 'No transpose', 'No transpose', I2, J2, + CALL SGEMM( 'No transpose', 'No transpose', I2, + $ J2, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+1, J+JB ), LDAB-1 ) @@ -411,7 +425,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A32 * - CALL SGEMM( 'No transpose', 'No transpose', I3, J2, + CALL SGEMM( 'No transpose', 'No transpose', I3, + $ J2, $ JB, -ONE, WORK31, LDWORK, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) @@ -431,7 +446,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A13 in the work array * - CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL STRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, $ WORK13, LDWORK ) * @@ -439,7 +455,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A23 * - CALL SGEMM( 'No transpose', 'No transpose', I2, J3, + CALL SGEMM( 'No transpose', 'No transpose', I2, + $ J3, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), $ LDAB-1 ) @@ -449,7 +466,8 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A33 * - CALL SGEMM( 'No transpose', 'No transpose', I3, J3, + CALL SGEMM( 'No transpose', 'No transpose', I3, + $ J3, $ JB, -ONE, WORK31, LDWORK, WORK13, $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF diff --git a/lapack-netlib/SRC/zgbtf2.c b/lapack-netlib/SRC/zgbtf2.c index 8d20eda268..0c0818541a 100644 --- a/lapack-netlib/SRC/zgbtf2.c +++ b/lapack-netlib/SRC/zgbtf2.c @@ -258,8 +258,6 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static doublecomplex c_b1 = {1.,0.}; @@ -273,7 +271,6 @@ e algorithm. */ /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download ZGBTF2 + dependencies */ /* > */ @@ -284,7 +281,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -382,9 +378,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complex16GBcomputational */ +/* > \ingroup gbtf2 */ /* > \par Further Details: */ /* ===================== */ @@ -418,21 +412,21 @@ f"> */ doublecomplex z__1; /* Local variables */ - integer i__, j; + integer i__, j, km, jp, ju, kv; + doublereal sfmin; extern /* Subroutine */ void zscal_(integer *, doublecomplex *, doublecomplex *, integer *), zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *), zswap_(integer *, doublecomplex *, integer *, doublecomplex *, integer *); - integer km, jp, ju, kv; + extern doublereal dlamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer izamax_(integer *, doublecomplex *, 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 */ /* ===================================================================== */ @@ -443,7 +437,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -476,6 +470,10 @@ f"> */ return; } +/* Compute machine safe minimum */ + + sfmin = dlamch_("S"); + /* Gaussian elimination with partial pivoting */ /* Set fill-in elements in columns KU+2 to KV to zero. */ @@ -540,14 +538,25 @@ f"> */ /* Compute multipliers. */ - z_div(&z__1, &c_b1, &ab[kv + 1 + j * ab_dim1]); - zscal_(&km, &z__1, &ab[kv + 2 + j * ab_dim1], &c__1); + if (z_abs(&ab[kv + 1 + j * ab_dim1]) >= sfmin) { + z_div(&z__1, &c_b1, &ab[kv + 1 + j * ab_dim1]); + zscal_(&km, &z__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; + z_div(&z__1, &ab[kv + 1 + i__ + j * ab_dim1], &ab[kv + + 1 + j * ab_dim1]); + ab[i__3].r = z__1.r, ab[i__3].i = z__1.i; +/* L35: */ + } + } /* Update trailing submatrix within the band. */ if (ju > j) { i__2 = ju - j; - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__3 = *ldab - 1; i__4 = *ldab - 1; zgeru_(&km, &i__2, &z__1, &ab[kv + 2 + j * ab_dim1], & diff --git a/lapack-netlib/SRC/zgbtf2.f b/lapack-netlib/SRC/zgbtf2.f index fbbe632476..b8f200910c 100644 --- a/lapack-netlib/SRC/zgbtf2.f +++ b/lapack-netlib/SRC/zgbtf2.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZGBTF2 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16GBcomputational +*> \ingroup gbtf2 * *> \par Further Details: * ===================== @@ -142,6 +140,7 @@ *> * ===================================================================== SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -163,11 +162,13 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) $ ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Local Scalars .. + DOUBLE PRECISION SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IZAMAX - EXTERNAL IZAMAX + EXTERNAL DLAMCH, IZAMAX * .. * .. External Subroutines .. EXTERNAL XERBLA, ZGERU, ZSCAL, ZSWAP @@ -206,6 +207,10 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -249,7 +254,14 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL ZSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL ZSCAL( 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. * diff --git a/lapack-netlib/SRC/zgbtrf.c b/lapack-netlib/SRC/zgbtrf.c index c94cac4a9d..00fc9ef76a 100644 --- a/lapack-netlib/SRC/zgbtrf.c +++ b/lapack-netlib/SRC/zgbtrf.c @@ -258,8 +258,6 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static doublecomplex c_b1 = {1.,0.}; @@ -273,7 +271,6 @@ static integer c__65 = 65; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download ZGBTRF + dependencies */ /* > */ @@ -284,7 +281,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -382,9 +378,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complex16GBcomputational */ +/* > \ingroup gbtrf */ /* > \par Further Details: */ /* ===================== */ @@ -417,8 +411,10 @@ f"> */ doublecomplex z__1; /* Local variables */ + integer i__, j, i2, i3, j2, j3, k2, jb, nb, ii, jj, jm, ip, jp, km, ju, + kv, nw; doublecomplex temp; - integer i__, j; + doublereal sfmin; extern /* Subroutine */ void zscal_(integer *, doublecomplex *, doublecomplex *, integer *), zgemm_(char *, char *, integer *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, @@ -426,7 +422,6 @@ f"> */ integer *); doublecomplex work13[4160] /* was [65][64] */, work31[4160] /* was [65][64] */; - integer i2, i3, j2, j3, k2; extern /* Subroutine */ void zgeru_(integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *), zcopy_(integer *, doublecomplex *, @@ -437,7 +432,7 @@ f"> */ integer *), zgbtf2_(integer *, integer *, integer *, integer *, doublecomplex *, integer *, integer *, integer *); - integer jb, nb, ii, jj, jm, ip, jp, km, ju, kv, nw; + extern doublereal dlamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen), izamax_(integer *, @@ -446,10 +441,9 @@ f"> */ 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 */ /* ===================================================================== */ @@ -460,7 +454,7 @@ f"> */ /* Parameter adjustments */ ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; + ab_offset = 1 + ab_dim1; ab -= ab_offset; --ipiv; @@ -493,6 +487,10 @@ f"> */ return; } +/* Compute machine safe minimum */ + + sfmin = dlamch_("S"); + /* Determine the block size for this environment */ nb = ilaenv_(&c__1, "ZGBTRF", " ", m, n, kl, ku, (ftnlen)6, (ftnlen)1); @@ -648,8 +646,19 @@ f"> */ /* Compute multipliers */ - z_div(&z__1, &c_b1, &ab[kv + 1 + jj * ab_dim1]); - zscal_(&km, &z__1, &ab[kv + 2 + jj * ab_dim1], &c__1); + if (z_abs(&ab[kv + 1 + jj * ab_dim1]) >= sfmin) { + z_div(&z__1, &c_b1, &ab[kv + 1 + jj * ab_dim1]); + zscal_(&km, &z__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; + z_div(&z__1, &ab[kv + 1 + i__ + jj * ab_dim1], & + ab[kv + 1 + jj * ab_dim1]); + ab[i__5].r = z__1.r, ab[i__5].i = z__1.i; +/* L75: */ + } + } /* Update trailing submatrix within the band and within */ /* the current block. JM is the index of the last column */ @@ -660,7 +669,7 @@ f"> */ jm = f2cmin(i__4,i__5); if (jm > jj) { i__4 = jm - jj; - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__5 = *ldab - 1; i__6 = *ldab - 1; zgeru_(&km, &i__4, &z__1, &ab[kv + 2 + jj * ab_dim1], @@ -754,7 +763,7 @@ f"> */ /* Update A22 */ - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__3 = *ldab - 1; i__4 = *ldab - 1; i__5 = *ldab - 1; @@ -769,7 +778,7 @@ f"> */ /* Update A32 */ - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__3 = *ldab - 1; i__4 = *ldab - 1; zgemm_("No transpose", "No transpose", &i3, &j2, &jb, @@ -808,7 +817,7 @@ f"> */ /* Update A23 */ - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__3 = *ldab - 1; i__4 = *ldab - 1; zgemm_("No transpose", "No transpose", &i2, &j3, &jb, @@ -821,7 +830,7 @@ f"> */ /* Update A33 */ - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; i__3 = *ldab - 1; zgemm_("No transpose", "No transpose", &i3, &j3, &jb, &z__1, work31, &c__65, work13, &c__65, &c_b1, diff --git a/lapack-netlib/SRC/zgbtrf.f b/lapack-netlib/SRC/zgbtrf.f index 83139b0919..9c8f55f7f3 100644 --- a/lapack-netlib/SRC/zgbtrf.f +++ b/lapack-netlib/SRC/zgbtrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZGBTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -115,7 +113,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16GBcomputational +*> \ingroup gbtrf * *> \par Further Details: * ===================== @@ -141,6 +139,7 @@ *> * ===================================================================== SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -167,17 +166,20 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW COMPLEX*16 TEMP + DOUBLE PRECISION SFMIN * .. * .. Local Arrays .. COMPLEX*16 WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER ILAENV, IZAMAX - EXTERNAL ILAENV, IZAMAX + EXTERNAL DLAMCH, ILAENV, IZAMAX * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZCOPY, ZGBTF2, ZGEMM, ZGERU, ZLASWP, + EXTERNAL XERBLA, ZCOPY, ZGBTF2, ZGEMM, ZGERU, + $ ZLASWP, $ ZSCAL, ZSWAP, ZTRSM * .. * .. Intrinsic Functions .. @@ -214,6 +216,10 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'ZGBTRF', ' ', M, N, KL, KU ) @@ -325,8 +331,15 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL ZSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL ZSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column @@ -394,7 +407,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A12 * - CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL ZTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J2, ONE, AB( KV+1, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1 ) * @@ -402,7 +416,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A22 * - CALL ZGEMM( 'No transpose', 'No transpose', I2, J2, + CALL ZGEMM( 'No transpose', 'No transpose', I2, + $ J2, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+1, J+JB ), LDAB-1 ) @@ -412,7 +427,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A32 * - CALL ZGEMM( 'No transpose', 'No transpose', I3, J2, + CALL ZGEMM( 'No transpose', 'No transpose', I3, + $ J2, $ JB, -ONE, WORK31, LDWORK, $ AB( KV+1-JB, J+JB ), LDAB-1, ONE, $ AB( KV+KL+1-JB, J+JB ), LDAB-1 ) @@ -432,7 +448,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A13 in the work array * - CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', + CALL ZTRSM( 'Left', 'Lower', 'No transpose', + $ 'Unit', $ JB, J3, ONE, AB( KV+1, J ), LDAB-1, $ WORK13, LDWORK ) * @@ -440,7 +457,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A23 * - CALL ZGEMM( 'No transpose', 'No transpose', I2, J3, + CALL ZGEMM( 'No transpose', 'No transpose', I2, + $ J3, $ JB, -ONE, AB( KV+1+JB, J ), LDAB-1, $ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ), $ LDAB-1 ) @@ -450,7 +468,8 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Update A33 * - CALL ZGEMM( 'No transpose', 'No transpose', I3, J3, + CALL ZGEMM( 'No transpose', 'No transpose', I3, + $ J3, $ JB, -ONE, WORK31, LDWORK, WORK13, $ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 ) END IF diff --git a/lapack-netlib/TESTING/LIN/alahd.f b/lapack-netlib/TESTING/LIN/alahd.f index 880effbebd..3f3930cdf1 100644 --- a/lapack-netlib/TESTING/LIN/alahd.f +++ b/lapack-netlib/TESTING/LIN/alahd.f @@ -831,7 +831,8 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) $ '2. First column zero', 15X, '6. Random, CNDNUM = .01/EPS', $ / 4X, '3. Last column zero', 16X, $ '7. Scaled near underflow', / 4X, - $ '4. Last n/2 columns zero', 11X, '8. Scaled near overflow' ) + $ '4. Last n/2 columns zero', 11X, '8. Scaled near overflow', + $ / 39X, '9. Scaled into the subnormal range' ) * * GT matrix types * diff --git a/lapack-netlib/TESTING/LIN/cchkaa.F b/lapack-netlib/TESTING/LIN/cchkaa.F index fa40000506..5e910c2632 100644 --- a/lapack-netlib/TESTING/LIN/cchkaa.F +++ b/lapack-netlib/TESTING/LIN/cchkaa.F @@ -204,19 +204,40 @@ PROGRAM CCHKAA * .. Allocate memory dynamically .. * ALLOCATE ( IWORK( 34*NMAX ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( A( ( KDMAX+1 )*NMAX, 8 ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( B( NMAX*MAXRHS, 4 ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( WORK( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( E( NMAX ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF * .. * .. Executable Statements .. * @@ -400,7 +421,7 @@ PROGRAM CCHKAA * IF( FATAL ) THEN WRITE( NOUT, FMT = 9999 ) - STOP + STOP 1 END IF * * Calculate and print the machine dependent constants. @@ -490,7 +511,7 @@ PROGRAM CCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1267,7 +1288,8 @@ PROGRAM CCHKAA DEALLOCATE (E, STAT = AllocateStatus) DEALLOCATE (S, STAT = AllocateStatus) DEALLOCATE (WORK, STAT = AllocateStatus) - DEALLOCATE (RWORK, STAT = AllocateStatus) + DEALLOCATE (RWORK, STAT = AllocateStatus) + DEALLOCATE (IWORK, STAT = AllocateStatus) * 9999 FORMAT( / ' Execution not attempted due to input errors' ) 9998 FORMAT( / ' End of tests' ) diff --git a/lapack-netlib/TESTING/LIN/cchkgb.f b/lapack-netlib/TESTING/LIN/cchkgb.f index 9bb6f29b79..c93abfa0ee 100644 --- a/lapack-netlib/TESTING/LIN/cchkgb.f +++ b/lapack-netlib/TESTING/LIN/cchkgb.f @@ -188,6 +188,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B, $ X, XACT, WORK, RWORK, IWORK, NOUT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -213,7 +214,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -227,6 +228,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + REAL ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -235,14 +237,17 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH REAL CLANGB, CLANGE, SGET06 EXTERNAL CLANGB, CLANGE, SGET06 + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CCOPY, CERRGE, CGBCON, $ CGBRFS, CGBT01, CGBT02, CGBT05, CGBTRF, CGBTRS, $ CGET04, CLACPY, CLARHS, CLASET, CLATB4, CLATMS, - $ XLAENV + $ CSSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC CMPLX, MAX, MIN @@ -400,6 +405,17 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. CLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = SLAMCH( 'Safe minimum' ) / 8 + CALL CSSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -477,13 +493,31 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL CGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = CLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( SISNAN( ANRMF ) .OR. + $ ANRMF.GT.SLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL CGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ SISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -495,7 +529,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = CLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/lapack-netlib/TESTING/LIN/dchkaa.F b/lapack-netlib/TESTING/LIN/dchkaa.F index ed5c7ab908..dd5110a820 100644 --- a/lapack-netlib/TESTING/LIN/dchkaa.F +++ b/lapack-netlib/TESTING/LIN/dchkaa.F @@ -507,7 +507,7 @@ PROGRAM DCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/lapack-netlib/TESTING/LIN/dchkgb.f b/lapack-netlib/TESTING/LIN/dchkgb.f index fa8dce5cff..7bb238a6c7 100644 --- a/lapack-netlib/TESTING/LIN/dchkgb.f +++ b/lapack-netlib/TESTING/LIN/dchkgb.f @@ -188,6 +188,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B, $ X, XACT, WORK, RWORK, IWORK, NOUT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -212,7 +213,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -226,6 +227,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN DOUBLE PRECISION AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + DOUBLE PRECISION ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -234,14 +236,17 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH DOUBLE PRECISION DGET06, DLANGB, DLANGE EXTERNAL DGET06, DLANGB, DLANGE + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DERRGE, DGBCON, $ DGBRFS, DGBT01, DGBT02, DGBT05, DGBTRF, DGBTRS, $ DGET04, DLACPY, DLARHS, DLASET, DLATB4, DLATMS, - $ XLAENV + $ DSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -400,6 +405,17 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. DLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = DLAMCH( 'Safe minimum' ) / 8 + CALL DSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -477,13 +493,31 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = DLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( DISNAN( ANRMF ) .OR. + $ ANRMF.GT.DLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL DGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ DISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -495,7 +529,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = DLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/lapack-netlib/TESTING/LIN/schkaa.F b/lapack-netlib/TESTING/LIN/schkaa.F index bb5b0af2cc..22694325d1 100644 --- a/lapack-netlib/TESTING/LIN/schkaa.F +++ b/lapack-netlib/TESTING/LIN/schkaa.F @@ -503,7 +503,7 @@ PROGRAM SCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/lapack-netlib/TESTING/LIN/schkgb.f b/lapack-netlib/TESTING/LIN/schkgb.f index 9ab8c0e018..b02ff97d26 100644 --- a/lapack-netlib/TESTING/LIN/schkgb.f +++ b/lapack-netlib/TESTING/LIN/schkgb.f @@ -188,6 +188,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B, $ X, XACT, WORK, RWORK, IWORK, NOUT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -212,7 +213,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -226,6 +227,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + REAL ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -234,14 +236,17 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH REAL SGET06, SLANGB, SLANGE EXTERNAL SGET06, SLANGB, SLANGE + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRGE, SGBCON, $ SGBRFS, SGBT01, SGBT02, SGBT05, SGBTRF, SGBTRS, $ SGET04, SLACPY, SLARHS, SLASET, SLATB4, SLATMS, - $ XLAENV + $ SSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -400,6 +405,17 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. SLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = SLAMCH( 'Safe minimum' ) / 8 + CALL SSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -477,13 +493,31 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = SLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( SISNAN( ANRMF ) .OR. + $ ANRMF.GT.SLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ SISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -495,7 +529,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/lapack-netlib/TESTING/LIN/zchkaa.F b/lapack-netlib/TESTING/LIN/zchkaa.F index a86793907d..50ab1103f7 100644 --- a/lapack-netlib/TESTING/LIN/zchkaa.F +++ b/lapack-netlib/TESTING/LIN/zchkaa.F @@ -204,19 +204,40 @@ PROGRAM ZCHKAA * .. Allocate memory dynamically .. * ALLOCATE ( IWORK( 34*NMAX ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( A ( (KDMAX+1) * NMAX, 8 ), STAT = AllocateStatus) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( B ( NMAX * MAXRHS, 4 ), STAT = AllocateStatus) - IF (AllocateStatus /= 0 ) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0 ) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( WORK ( NMAX, NMAX+MAXRHS+10 ), STAT = AllocateStatus) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( E( NMAX ), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( S( 2*NMAX ), STAT = AllocateStatus) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( RWORK( 150*NMAX+2*MAXRHS ), STAT = AllocateStatus) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF * .. * .. Executable Statements .. * @@ -400,7 +421,7 @@ PROGRAM ZCHKAA * IF( FATAL ) THEN WRITE( NOUT, FMT = 9999 ) - STOP + STOP 1 END IF * * Calculate and print the machine dependent constants. @@ -490,7 +511,7 @@ PROGRAM ZCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1300,8 +1321,9 @@ PROGRAM ZCHKAA DEALLOCATE (B, STAT = AllocateStatus) DEALLOCATE (E, STAT = AllocateStatus) DEALLOCATE (S, STAT = AllocateStatus) + DEALLOCATE (WORK, STAT = AllocateStatus) DEALLOCATE (RWORK, STAT = AllocateStatus) - DEALLOCATE (WORK, STAT = AllocateStatus) + DEALLOCATE (IWORK, STAT = AllocateStatus) * 9999 FORMAT( / ' Execution not attempted due to input errors' ) 9998 FORMAT( / ' End of tests' ) diff --git a/lapack-netlib/TESTING/LIN/zchkgb.f b/lapack-netlib/TESTING/LIN/zchkgb.f index 3b0cd707ab..b4a2883626 100644 --- a/lapack-netlib/TESTING/LIN/zchkgb.f +++ b/lapack-netlib/TESTING/LIN/zchkgb.f @@ -188,6 +188,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NSVAL, THRESH, TSTERR, A, LA, AFAC, LAFAC, B, $ X, XACT, WORK, RWORK, IWORK, NOUT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -213,7 +214,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -227,6 +228,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN DOUBLE PRECISION AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + DOUBLE PRECISION ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -235,14 +237,17 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH DOUBLE PRECISION DGET06, ZLANGB, ZLANGE EXTERNAL DGET06, ZLANGB, ZLANGE + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. - EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZCOPY, ZERRGE, - $ ZGBCON, ZGBRFS, ZGBT01, ZGBT02, ZGBT05, ZGBTRF, - $ ZGBTRS, ZGET04, ZLACPY, ZLARHS, ZLASET, ZLATB4, - $ ZLATMS + EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZCOPY, ZDSCAL, + $ ZERRGE, ZGBCON, ZGBRFS, ZGBT01, ZGBT02, ZGBT05, + $ ZGBTRF, ZGBTRS, ZGET04, ZLACPY, ZLARHS, ZLASET, + $ ZLATB4, ZLATMS * .. * .. Intrinsic Functions .. INTRINSIC DCMPLX, MAX, MIN @@ -400,6 +405,17 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. ZLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = DLAMCH( 'Safe minimum' ) / 8 + CALL ZDSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -477,13 +493,31 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = ZLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( DISNAN( ANRMF ) .OR. + $ ANRMF.GT.DLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL ZGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ DISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -495,7 +529,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = ZLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/lapack-netlib/TESTING/ctest.in b/lapack-netlib/TESTING/ctest.in index 4e30224d73..d740caddfa 100644 --- a/lapack-netlib/TESTING/ctest.in +++ b/lapack-netlib/TESTING/ctest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits CGE 11 List types on next line if 0 < NTYPES < 11 -CGB 8 List types on next line if 0 < NTYPES < 8 +CGB 9 List types on next line if 0 < NTYPES < 9 CGT 12 List types on next line if 0 < NTYPES < 12 CPO 9 List types on next line if 0 < NTYPES < 9 CPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/lapack-netlib/TESTING/dtest.in b/lapack-netlib/TESTING/dtest.in index 34a6dd14b2..71ef5be566 100644 --- a/lapack-netlib/TESTING/dtest.in +++ b/lapack-netlib/TESTING/dtest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits DGE 11 List types on next line if 0 < NTYPES < 11 -DGB 8 List types on next line if 0 < NTYPES < 8 +DGB 9 List types on next line if 0 < NTYPES < 9 DGT 12 List types on next line if 0 < NTYPES < 12 DPO 9 List types on next line if 0 < NTYPES < 9 DPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/lapack-netlib/TESTING/stest.in b/lapack-netlib/TESTING/stest.in index b95b4ed322..839ca2ee9b 100644 --- a/lapack-netlib/TESTING/stest.in +++ b/lapack-netlib/TESTING/stest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits SGE 11 List types on next line if 0 < NTYPES < 11 -SGB 8 List types on next line if 0 < NTYPES < 8 +SGB 9 List types on next line if 0 < NTYPES < 9 SGT 12 List types on next line if 0 < NTYPES < 12 SPO 9 List types on next line if 0 < NTYPES < 9 SPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/lapack-netlib/TESTING/ztest.in b/lapack-netlib/TESTING/ztest.in index bf4c9d100b..eeb1258b7f 100644 --- a/lapack-netlib/TESTING/ztest.in +++ b/lapack-netlib/TESTING/ztest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits ZGE 11 List types on next line if 0 < NTYPES < 11 -ZGB 8 List types on next line if 0 < NTYPES < 8 +ZGB 9 List types on next line if 0 < NTYPES < 9 ZGT 12 List types on next line if 0 < NTYPES < 12 ZPO 9 List types on next line if 0 < NTYPES < 9 ZPS 9 List types on next line if 0 < NTYPES < 9