Skip to content

Commit b8e55b1

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #6175 from martin-frbg/lapack1377
Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite factorizations (Reference-LAPACK PR 1377)
2 parents 12f1ec9 + 6b8e4e5 commit b8e55b1

67 files changed

Lines changed: 2448 additions & 4183 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎lapack-netlib/SRC/chetf2.c‎

Lines changed: 43 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -249,7 +249,6 @@ typedef struct Namelist Namelist;
249249

250250

251251

252-
253252
/* Table of constant values */
254253

255254
static integer c__1 = 1;
@@ -262,7 +261,6 @@ thod (unblocked algorithm calling Level 2 BLAS). */
262261
/* Online html documentation available at */
263262
/* http://www.netlib.org/lapack/explore-html/ */
264263

265-
/* > \htmlonly */
266264
/* > Download CHETF2 + dependencies */
267265
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2.
268266
f"> */
@@ -273,7 +271,6 @@ f"> */
273271
/* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2.
274272
f"> */
275273
/* > [TXT]</a> */
276-
/* > \endhtmlonly */
277274

278275
/* Definition: */
279276
/* =========== */
@@ -383,9 +380,7 @@ f"> */
383380
/* > \author Univ. of Colorado Denver */
384381
/* > \author NAG Ltd. */
385382

386-
/* > \date December 2016 */
387-
388-
/* > \ingroup complexHEcomputational */
383+
/* > \ingroup hetf2 */
389384

390385
/* > \par Further Details: */
391386
/* ===================== */
@@ -445,45 +440,43 @@ f"> */
445440
{
446441
/* System generated locals */
447442
integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
448-
real r__1=0., r__2=0., r__3=0., r__4=0.;
449-
complex q__1={0.,0.}, q__2={0.,0.}, q__3={0.,0.}, q__4={0.,0.}, q__5={0.,0.}, q__6={0.,0.};
443+
real r__1, r__2, r__3, r__4;
444+
complex q__1, q__2, q__3, q__4, q__5, q__6;
450445

451446
/* Local variables */
452-
extern /* Subroutine */ void cher_(char *, integer *, real *, complex *,
453-
integer *, complex *, integer *);
454-
integer imax, jmax;
455447
real d__;
456448
integer i__, j, k;
457449
complex t;
450+
real r1, d11;
451+
complex d12;
452+
real d22;
453+
complex d21;
454+
integer kk, kp;
455+
complex wk;
456+
real tt;
457+
complex wkm1, wkp1;
458+
extern /* Subroutine */ void cher_(char *, integer *, real *, complex *,
459+
integer *, complex *, integer *);
460+
integer imax, jmax;
458461
real alpha;
459462
extern logical lsame_(char *, char *);
460463
extern /* Subroutine */ void cswap_(integer *, complex *, integer *,
461464
complex *, integer *);
462465
integer kstep;
463466
logical upper;
464-
real r1, d11;
465-
complex d12;
466-
real d22;
467-
complex d21;
468467
extern real slapy2_(real *, real *);
469-
integer kk, kp;
470468
real absakk;
471-
complex wk;
472469
extern integer icamax_(integer *, complex *, integer *);
473-
real tt;
474470
extern /* Subroutine */ void csscal_(integer *, real *, complex *, integer
475-
*);
476-
extern void xerbla_(char *, integer *, ftnlen);
471+
*), xerbla_(char *, integer *, ftnlen);
477472
real colmax;
478473
extern logical sisnan_(real *);
479474
real rowmax;
480-
complex wkm1, wkp1;
481475

482476

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

488481

489482
/* ===================================================================== */
@@ -493,7 +486,7 @@ f"> */
493486

494487
/* Parameter adjustments */
495488
a_dim1 = *lda;
496-
a_offset = 1 + a_dim1 * 1;
489+
a_offset = 1 + a_dim1;
497490
a -= a_offset;
498491
--ipiv;
499492

@@ -729,27 +722,28 @@ f"> */
729722
i__1 = k - 1 + k * a_dim1;
730723
q__1.r = a[i__1].r / d__, q__1.i = a[i__1].i / d__;
731724
d12.r = q__1.r, d12.i = q__1.i;
732-
d__ = tt / d__;
733725

734726
for (j = k - 2; j >= 1; --j) {
735727
i__1 = j + (k - 1) * a_dim1;
736-
q__3.r = d11 * a[i__1].r, q__3.i = d11 * a[i__1].i;
737-
r_cnjg(&q__5, &d12);
728+
q__4.r = d11 * a[i__1].r, q__4.i = d11 * a[i__1].i;
729+
r_cnjg(&q__6, &d12);
738730
i__2 = j + k * a_dim1;
739-
q__4.r = q__5.r * a[i__2].r - q__5.i * a[i__2].i,
740-
q__4.i = q__5.r * a[i__2].i + q__5.i * a[i__2]
731+
q__5.r = q__6.r * a[i__2].r - q__6.i * a[i__2].i,
732+
q__5.i = q__6.r * a[i__2].i + q__6.i * a[i__2]
741733
.r;
742-
q__2.r = q__3.r - q__4.r, q__2.i = q__3.i - q__4.i;
743-
q__1.r = d__ * q__2.r, q__1.i = d__ * q__2.i;
734+
q__3.r = q__4.r - q__5.r, q__3.i = q__4.i - q__5.i;
735+
q__2.r = q__3.r / d__, q__2.i = q__3.i / d__;
736+
q__1.r = tt * q__2.r, q__1.i = tt * q__2.i;
744737
wkm1.r = q__1.r, wkm1.i = q__1.i;
745738
i__1 = j + k * a_dim1;
746-
q__3.r = d22 * a[i__1].r, q__3.i = d22 * a[i__1].i;
739+
q__4.r = d22 * a[i__1].r, q__4.i = d22 * a[i__1].i;
747740
i__2 = j + (k - 1) * a_dim1;
748-
q__4.r = d12.r * a[i__2].r - d12.i * a[i__2].i,
749-
q__4.i = d12.r * a[i__2].i + d12.i * a[i__2]
741+
q__5.r = d12.r * a[i__2].r - d12.i * a[i__2].i,
742+
q__5.i = d12.r * a[i__2].i + d12.i * a[i__2]
750743
.r;
751-
q__2.r = q__3.r - q__4.r, q__2.i = q__3.i - q__4.i;
752-
q__1.r = d__ * q__2.r, q__1.i = d__ * q__2.i;
744+
q__3.r = q__4.r - q__5.r, q__3.i = q__4.i - q__5.i;
745+
q__2.r = q__3.r / d__, q__2.i = q__3.i / d__;
746+
q__1.r = tt * q__2.r, q__1.i = tt * q__2.i;
753747
wk.r = q__1.r, wk.i = q__1.i;
754748
for (i__ = j; i__ >= 1; --i__) {
755749
i__1 = i__ + j * a_dim1;
@@ -1017,28 +1011,29 @@ f"> */
10171011
i__1 = k + 1 + k * a_dim1;
10181012
q__1.r = a[i__1].r / d__, q__1.i = a[i__1].i / d__;
10191013
d21.r = q__1.r, d21.i = q__1.i;
1020-
d__ = tt / d__;
10211014

10221015
i__1 = *n;
10231016
for (j = k + 2; j <= i__1; ++j) {
10241017
i__2 = j + k * a_dim1;
1025-
q__3.r = d11 * a[i__2].r, q__3.i = d11 * a[i__2].i;
1018+
q__4.r = d11 * a[i__2].r, q__4.i = d11 * a[i__2].i;
10261019
i__3 = j + (k + 1) * a_dim1;
1027-
q__4.r = d21.r * a[i__3].r - d21.i * a[i__3].i,
1028-
q__4.i = d21.r * a[i__3].i + d21.i * a[i__3]
1020+
q__5.r = d21.r * a[i__3].r - d21.i * a[i__3].i,
1021+
q__5.i = d21.r * a[i__3].i + d21.i * a[i__3]
10291022
.r;
1030-
q__2.r = q__3.r - q__4.r, q__2.i = q__3.i - q__4.i;
1031-
q__1.r = d__ * q__2.r, q__1.i = d__ * q__2.i;
1023+
q__3.r = q__4.r - q__5.r, q__3.i = q__4.i - q__5.i;
1024+
q__2.r = q__3.r / d__, q__2.i = q__3.i / d__;
1025+
q__1.r = tt * q__2.r, q__1.i = tt * q__2.i;
10321026
wk.r = q__1.r, wk.i = q__1.i;
10331027
i__2 = j + (k + 1) * a_dim1;
1034-
q__3.r = d22 * a[i__2].r, q__3.i = d22 * a[i__2].i;
1035-
r_cnjg(&q__5, &d21);
1028+
q__4.r = d22 * a[i__2].r, q__4.i = d22 * a[i__2].i;
1029+
r_cnjg(&q__6, &d21);
10361030
i__3 = j + k * a_dim1;
1037-
q__4.r = q__5.r * a[i__3].r - q__5.i * a[i__3].i,
1038-
q__4.i = q__5.r * a[i__3].i + q__5.i * a[i__3]
1031+
q__5.r = q__6.r * a[i__3].r - q__6.i * a[i__3].i,
1032+
q__5.i = q__6.r * a[i__3].i + q__6.i * a[i__3]
10391033
.r;
1040-
q__2.r = q__3.r - q__4.r, q__2.i = q__3.i - q__4.i;
1041-
q__1.r = d__ * q__2.r, q__1.i = d__ * q__2.i;
1034+
q__3.r = q__4.r - q__5.r, q__3.i = q__4.i - q__5.i;
1035+
q__2.r = q__3.r / d__, q__2.i = q__3.i / d__;
1036+
q__1.r = tt * q__2.r, q__1.i = tt * q__2.i;
10421037
wkp1.r = q__1.r, wkp1.i = q__1.i;
10431038
i__2 = *n;
10441039
for (i__ = j; i__ <= i__2; ++i__) {

‎lapack-netlib/SRC/chetf2.f‎

Lines changed: 16 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download CHETF2 + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetf2.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -127,7 +125,7 @@
127125
*> \author Univ. of Colorado Denver
128126
*> \author NAG Ltd.
129127
*
130-
*> \ingroup complexHEcomputational
128+
*> \ingroup hetf2
131129
*
132130
*> \par Further Details:
133131
* =====================
@@ -183,6 +181,7 @@
183181
*>
184182
* =====================================================================
185183
SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
184+
IMPLICIT NONE
186185
*
187186
* -- LAPACK computational routine --
188187
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -284,7 +283,8 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
284283
COLMAX = ZERO
285284
END IF
286285
*
287-
IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
286+
IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR.
287+
$ SISNAN(ABSAKK) ) THEN
288288
*
289289
* Column K is or underflow, or contains a NaN:
290290
* set INFO and continue
@@ -403,11 +403,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
403403
D11 = REAL( A( K, K ) ) / D
404404
TT = ONE / ( D11*D22-ONE )
405405
D12 = A( K-1, K ) / D
406-
D = TT / D
407406
*
408407
DO 40 J = K - 2, 1, -1
409-
WKM1 = D*( D11*A( J, K-1 )-CONJG( D12 )*A( J, K ) )
410-
WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
408+
WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )*
409+
$ A( J, K ) ) / D )
410+
WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D )
411411
DO 30 I = J, 1, -1
412412
A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
413413
$ A( I, K-1 )*CONJG( WKM1 )
@@ -468,7 +468,8 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
468468
COLMAX = ZERO
469469
END IF
470470
*
471-
IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
471+
IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR.
472+
$ SISNAN(ABSAKK) ) THEN
472473
*
473474
* Column K is zero or underflow, contains a NaN:
474475
* set INFO and continue
@@ -491,7 +492,8 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
491492
JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
492493
ROWMAX = CABS1( A( IMAX, JMAX ) )
493494
IF( IMAX.LT.N ) THEN
494-
JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
495+
JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ),
496+
$ 1 )
495497
ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
496498
END IF
497499
*
@@ -524,7 +526,8 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
524526
* submatrix A(k:n,k:n)
525527
*
526528
IF( KP.LT.N )
527-
$ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
529+
$ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ),
530+
$ 1 )
528531
DO 60 J = KK + 1, KP - 1
529532
T = CONJG( A( J, KK ) )
530533
A( J, KK ) = CONJG( A( KP, J ) )
@@ -590,11 +593,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
590593
D22 = REAL( A( K, K ) ) / D
591594
TT = ONE / ( D11*D22-ONE )
592595
D21 = A( K+1, K ) / D
593-
D = TT / D
594596
*
595597
DO 80 J = K + 2, N
596-
WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
597-
WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) )
598+
WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D )
599+
WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )*
600+
$ A( J, K ) ) / D )
598601
DO 70 I = J, N
599602
A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
600603
$ A( I, K+1 )*CONJG( WKP1 )

0 commit comments

Comments
 (0)