diff --git a/SRC/chetf2.f b/SRC/chetf2.f index 0a6adfcef..28dcbd3c3 100644 --- a/SRC/chetf2.f +++ b/SRC/chetf2.f @@ -403,11 +403,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = REAL( A( K, K ) ) / D TT = ONE / ( D11*D22-ONE ) D12 = A( K-1, K ) / D - D = TT / D * DO 40 J = K - 2, 1, -1 - WKM1 = D*( D11*A( J, K-1 )-CONJG( D12 )*A( J, K ) ) - WK = D*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) DO 30 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) - $ A( I, K-1 )*CONJG( WKM1 ) @@ -593,11 +593,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = REAL( A( K, K ) ) / D TT = ONE / ( D11*D22-ONE ) D21 = A( K+1, K ) / D - D = TT / D * DO 80 J = K + 2, N - WK = D*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )* + $ A( J, K ) ) / D ) DO 70 I = J, N A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) - $ A( I, K+1 )*CONJG( WKP1 ) diff --git a/SRC/chetf2_rk.f b/SRC/chetf2_rk.f index e9b427434..bcd8e1d0b 100644 --- a/SRC/chetf2_rk.f +++ b/SRC/chetf2_rk.f @@ -617,24 +617,24 @@ SUBROUTINE CHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )* - $ A( J, K ) ) - WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) * * Perform a rank-2 update of A(1:k-2,1:k-2) * DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*CONJG( WK ) - - $ ( A( I, K-1 ) / D )*CONJG( WKM1 ) + $ A( I, K )*CONJG( WK ) - + $ A( I, K-1 )*CONJG( WKM1 ) 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D - A( J, K-1 ) = WKM1 / D + A( J, K ) = WK + A( J, K-1 ) = WKM1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) * @@ -978,24 +978,24 @@ SUBROUTINE CHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )* - $ A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )* + $ A( J, K ) ) / D ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*CONJG( WK ) - - $ ( A( I, K+1 ) / D )*CONJG( WKP1 ) + $ A( I, K )*CONJG( WK ) - + $ A( I, K+1 )*CONJG( WKP1 ) 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D - A( J, K+1 ) = WKP1 / D + A( J, K ) = WK + A( J, K+1 ) = WKP1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) * diff --git a/SRC/chetf2_rook.f b/SRC/chetf2_rook.f index 1f49b604b..fcaef0922 100644 --- a/SRC/chetf2_rook.f +++ b/SRC/chetf2_rook.f @@ -536,24 +536,24 @@ SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )* - $ A( J, K ) ) - WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) * * Perform a rank-2 update of A(1:k-2,1:k-2) * DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*CONJG( WK ) - - $ ( A( I, K-1 ) / D )*CONJG( WKM1 ) + $ A( I, K )*CONJG( WK ) - + $ A( I, K-1 )*CONJG( WKM1 ) 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D - A( J, K-1 ) = WKM1 / D + A( J, K ) = WK + A( J, K-1 ) = WKM1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) * @@ -857,24 +857,24 @@ SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )* - $ A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )* + $ A( J, K ) ) / D ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*CONJG( WK ) - - $ ( A( I, K+1 ) / D )*CONJG( WKP1 ) + $ A( I, K )*CONJG( WK ) - + $ A( I, K+1 )*CONJG( WKP1 ) 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D - A( J, K+1 ) = WKP1 / D + A( J, K ) = WK + A( J, K+1 ) = WKP1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO ) * diff --git a/SRC/chptrf.f b/SRC/chptrf.f index 3d584fd51..bc4048f33 100644 --- a/SRC/chptrf.f +++ b/SRC/chptrf.f @@ -389,13 +389,12 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) D11 = REAL( AP( K+( K-1 )*K / 2 ) ) / D TT = ONE / ( D11*D22-ONE ) D12 = AP( K-1+( K-1 )*K / 2 ) / D - D = TT / D * DO 50 J = K - 2, 1, -1 - WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ CONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) - WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12* - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = TT*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ CONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) / D ) + WK = TT*( ( D22*AP( J+( K-1 )*K / 2 )-D12* + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*CONJG( WK ) - @@ -603,13 +602,13 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) D22 = REAL( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D TT = ONE / ( D11*D22-ONE ) D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D - D = TT / D * DO 100 J = K + 2, N - WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21* - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- - $ CONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 2 ) ) + WK = TT*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21* + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D ) + WKP1 = TT*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + $ CONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 2 ) + $ ) / D ) DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / diff --git a/SRC/clahef.f b/SRC/clahef.f index 21c9f0986..d1654b8b6 100644 --- a/SRC/clahef.f +++ b/SRC/clahef.f @@ -481,17 +481,17 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = * ( ( -1 ) ( D22 ) ) * -* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ), * ( ( -1 ) ( D22 ) ) * -* = ( conj(D21)*( D11 ) D21*( -1 ) ) -* ( ( -1 ) ( D22 ) ), -* * where D11 = d22/d21, * D22 = d11/conj(d21), -* D21 = T/d21, * T = 1/(D22*D11-1). * +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 or conj(d21) and then scaled by T. +* * (NOTE: No need to check for division by ZERO, * since that was ensured earlier in pivot search: * (a) d21 != 0, since in 2x2 pivot case(4) @@ -503,16 +503,16 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, D11 = W( K, KW ) / CONJG( D21 ) D22 = W( K-1, KW-1 ) / D21 T = ONE / ( REAL( D11*D22 )-ONE ) - D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = CONJG( D21 )* - $ ( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( ( D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ CONJG( D21 ) ) 20 CONTINUE END IF * @@ -828,17 +828,17 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = * ( ( -1 ) ( D22 ) ) * -* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ), * ( ( -1 ) ( D22 ) ) * -* = ( conj(D21)*( D11 ) D21*( -1 ) ) -* ( ( -1 ) ( D22 ) ) -* * where D11 = d22/d21, * D22 = d11/conj(d21), -* D21 = T/d21, * T = 1/(D22*D11-1). * +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 or conj(d21) and then scaled by T. +* * (NOTE: No need to check for division by ZERO, * since that was ensured earlier in pivot search: * (a) d21 != 0, since in 2x2 pivot case(4) @@ -850,16 +850,16 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / CONJG( D21 ) T = ONE / ( REAL( D11*D22 )-ONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = CONJG( D21 )* - $ ( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ CONJG( D21 ) ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/clasyf.f b/SRC/clasyf.f index 926c62de6..0b9248a29 100644 --- a/SRC/clasyf.f +++ b/SRC/clasyf.f @@ -432,8 +432,9 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K-1, KW ) D11 = W( K, KW ) / D21 @@ -444,10 +445,11 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * - D21 = T / D21 DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ D21 ) 20 CONTINUE END IF * @@ -712,22 +714,24 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K+1, K ) D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ D21 ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/csptrf.f b/SRC/csptrf.f index 64f8e678b..c7d2d7f41 100644 --- a/SRC/csptrf.f +++ b/SRC/csptrf.f @@ -375,13 +375,12 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 D11 = AP( K+( K-1 )*K / 2 ) / D12 T = CONE / ( D11*D22-CONE ) - D12 = T / D12 * DO 50 J = K - 2, 1, -1 - WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ AP( J+( K-1 )*K / 2 ) ) - WK = D12*( D22*AP( J+( K-1 )*K / 2 )- - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = T*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ AP( J+( K-1 )*K / 2 ) ) / D12 ) + WK = T*( ( D22*AP( J+( K-1 )*K / 2 )- + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D12 ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*WK - @@ -576,13 +575,12 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * DO 100 J = K + 2, N - WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- - $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) + WK = T*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D21 ) + WKP1 = T*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) / D21 ) DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / diff --git a/SRC/csytf2.f b/SRC/csytf2.f index 01c882fad..09228ec9f 100644 --- a/SRC/csytf2.f +++ b/SRC/csytf2.f @@ -394,11 +394,10 @@ SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = A( K-1, K-1 ) / D12 D11 = A( K, K ) / D12 T = CONE / ( D11*D22-CONE ) - D12 = T / D12 * DO 30 J = K - 2, 1, -1 - WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) - WK = D12*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K-1 )*WKM1 @@ -569,11 +568,10 @@ SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = A( K+1, K+1 ) / D21 D22 = A( K, K ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * DO 60 J = K + 2, N - WK = D21*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) DO 50 I = J, N A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K+1 )*WKP1 diff --git a/SRC/csytf2_rk.f b/SRC/csytf2_rk.f index a2ff8ad9e..a52f69a0b 100644 --- a/SRC/csytf2_rk.f +++ b/SRC/csytf2_rk.f @@ -582,18 +582,18 @@ SUBROUTINE CSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -898,22 +898,22 @@ SUBROUTINE CSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/csytf2_rook.f b/SRC/csytf2_rook.f index df504cef8..50f2894b4 100644 --- a/SRC/csytf2_rook.f +++ b/SRC/csytf2_rook.f @@ -502,18 +502,18 @@ SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -776,22 +776,22 @@ SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/dlasyf.f b/SRC/dlasyf.f index 68e0087f3..a3e66e965 100644 --- a/SRC/dlasyf.f +++ b/SRC/dlasyf.f @@ -424,22 +424,24 @@ SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K-1, KW ) D11 = W( K, KW ) / D21 D22 = W( K-1, KW-1 ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ D21 ) 20 CONTINUE END IF * @@ -703,22 +705,24 @@ SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K+1, K ) D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ D21 ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/dsptrf.f b/SRC/dsptrf.f index acd029043..2a05f0138 100644 --- a/SRC/dsptrf.f +++ b/SRC/dsptrf.f @@ -368,13 +368,12 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 D11 = AP( K+( K-1 )*K / 2 ) / D12 T = ONE / ( D11*D22-ONE ) - D12 = T / D12 * DO 50 J = K - 2, 1, -1 - WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ AP( J+( K-1 )*K / 2 ) ) - WK = D12*( D22*AP( J+( K-1 )*K / 2 )- - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = T*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ AP( J+( K-1 )*K / 2 ) ) / D12 ) + WK = T*( ( D22*AP( J+( K-1 )*K / 2 )- + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D12 ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*WK - @@ -570,13 +569,12 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * DO 100 J = K + 2, N - WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- - $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) + WK = T*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D21 ) + WKP1 = T*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) / D21 ) * DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* diff --git a/SRC/dsytf2.f b/SRC/dsytf2.f index bd4340a6a..39685c36a 100644 --- a/SRC/dsytf2.f +++ b/SRC/dsytf2.f @@ -390,11 +390,10 @@ SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = A( K-1, K-1 ) / D12 D11 = A( K, K ) / D12 T = ONE / ( D11*D22-ONE ) - D12 = T / D12 * DO 30 J = K - 2, 1, -1 - WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) - WK = D12*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K-1 )*WKM1 @@ -565,12 +564,11 @@ SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = A( K+1, K+1 ) / D21 D22 = A( K, K ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * DO 60 J = K + 2, N * - WK = D21*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * DO 50 I = J, N A( I, J ) = A( I, J ) - A( I, K )*WK - diff --git a/SRC/dsytf2_rk.f b/SRC/dsytf2_rk.f index 0c464ec0f..1628d9f0e 100644 --- a/SRC/dsytf2_rk.f +++ b/SRC/dsytf2_rk.f @@ -573,18 +573,18 @@ SUBROUTINE DSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -889,22 +889,22 @@ SUBROUTINE DSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/dsytf2_rook.f b/SRC/dsytf2_rook.f index 270deab60..e5935aec7 100644 --- a/SRC/dsytf2_rook.f +++ b/SRC/dsytf2_rook.f @@ -494,18 +494,18 @@ SUBROUTINE DSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -768,22 +768,22 @@ SUBROUTINE DSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/slasyf.f b/SRC/slasyf.f index 129607a2b..195a388d8 100644 --- a/SRC/slasyf.f +++ b/SRC/slasyf.f @@ -424,22 +424,24 @@ SUBROUTINE SLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K-1, KW ) D11 = W( K, KW ) / D21 D22 = W( K-1, KW-1 ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ D21 ) 20 CONTINUE END IF * @@ -703,22 +705,24 @@ SUBROUTINE SLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K+1, K ) D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ D21 ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/ssptrf.f b/SRC/ssptrf.f index 3c4456a14..dd54eb1f5 100644 --- a/SRC/ssptrf.f +++ b/SRC/ssptrf.f @@ -366,13 +366,12 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 D11 = AP( K+( K-1 )*K / 2 ) / D12 T = ONE / ( D11*D22-ONE ) - D12 = T / D12 * DO 50 J = K - 2, 1, -1 - WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ AP( J+( K-1 )*K / 2 ) ) - WK = D12*( D22*AP( J+( K-1 )*K / 2 )- - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = T*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ AP( J+( K-1 )*K / 2 ) ) / D12 ) + WK = T*( ( D22*AP( J+( K-1 )*K / 2 )- + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D12 ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*WK - @@ -568,13 +567,12 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * DO 100 J = K + 2, N - WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- - $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) + WK = T*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D21 ) + WKP1 = T*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) / D21 ) * DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* diff --git a/SRC/ssytf2.f b/SRC/ssytf2.f index d5defcccc..739852112 100644 --- a/SRC/ssytf2.f +++ b/SRC/ssytf2.f @@ -391,11 +391,10 @@ SUBROUTINE SSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = A( K-1, K-1 ) / D12 D11 = A( K, K ) / D12 T = ONE / ( D11*D22-ONE ) - D12 = T / D12 * DO 30 J = K - 2, 1, -1 - WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) - WK = D12*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K-1 )*WKM1 @@ -566,12 +565,11 @@ SUBROUTINE SSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = A( K+1, K+1 ) / D21 D22 = A( K, K ) / D21 T = ONE / ( D11*D22-ONE ) - D21 = T / D21 * DO 60 J = K + 2, N * - WK = D21*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * DO 50 I = J, N A( I, J ) = A( I, J ) - A( I, K )*WK - diff --git a/SRC/ssytf2_rk.f b/SRC/ssytf2_rk.f index d78f621f0..12e76adc8 100644 --- a/SRC/ssytf2_rk.f +++ b/SRC/ssytf2_rk.f @@ -573,18 +573,18 @@ SUBROUTINE SSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -889,22 +889,22 @@ SUBROUTINE SSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/ssytf2_rook.f b/SRC/ssytf2_rook.f index 9e6488728..dae4093a7 100644 --- a/SRC/ssytf2_rook.f +++ b/SRC/ssytf2_rook.f @@ -494,18 +494,18 @@ SUBROUTINE SSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -768,22 +768,22 @@ SUBROUTINE SSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/zhetf2.f b/SRC/zhetf2.f index eb47f4925..40a7b0398 100644 --- a/SRC/zhetf2.f +++ b/SRC/zhetf2.f @@ -418,12 +418,11 @@ SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = DBLE( A( K, K ) ) / D TT = ONE / ( D11*D22-ONE ) D12 = A( K-1, K ) / D - D = TT / D * DO 40 J = K - 2, 1, -1 - WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )* - $ A( J, K ) ) - WK = D*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-DCONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) DO 30 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) - $ A( I, K-1 )*DCONJG( WKM1 ) @@ -619,12 +618,11 @@ SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = DBLE( A( K, K ) ) / D TT = ONE / ( D11*D22-ONE ) D21 = A( K+1, K ) / D - D = TT / D * DO 80 J = K + 2, N - WK = D*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )* - $ A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-DCONJG( D21 )* + $ A( J, K ) ) / D ) DO 70 I = J, N A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) - $ A( I, K+1 )*DCONJG( WKP1 ) diff --git a/SRC/zhetf2_rk.f b/SRC/zhetf2_rk.f index 5e2afb066..b8c9bc905 100644 --- a/SRC/zhetf2_rk.f +++ b/SRC/zhetf2_rk.f @@ -617,24 +617,24 @@ SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WKM1 = TT*( D11*A( J, K-1 )-DCONJG( D12 )* - $ A( J, K ) ) - WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-DCONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) * * Perform a rank-2 update of A(1:k-2,1:k-2) * DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*DCONJG( WK ) - - $ ( A( I, K-1 ) / D )*DCONJG( WKM1 ) + $ A( I, K )*DCONJG( WK ) - + $ A( I, K-1 )*DCONJG( WKM1 ) 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D - A( J, K-1 ) = WKM1 / D + A( J, K ) = WK + A( J, K-1 ) = WKM1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO ) * @@ -978,24 +978,24 @@ SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = TT*( D22*A( J, K+1 )-DCONJG( D21 )* - $ A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-DCONJG( D21 )* + $ A( J, K ) ) / D ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*DCONJG( WK ) - - $ ( A( I, K+1 ) / D )*DCONJG( WKP1 ) + $ A( I, K )*DCONJG( WK ) - + $ A( I, K+1 )*DCONJG( WKP1 ) 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D - A( J, K+1 ) = WKP1 / D + A( J, K ) = WK + A( J, K+1 ) = WKP1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO ) * diff --git a/SRC/zhetf2_rook.f b/SRC/zhetf2_rook.f index 4e86aec7b..333150bcb 100644 --- a/SRC/zhetf2_rook.f +++ b/SRC/zhetf2_rook.f @@ -536,24 +536,24 @@ SUBROUTINE ZHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WKM1 = TT*( D11*A( J, K-1 )-DCONJG( D12 )* - $ A( J, K ) ) - WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) ) + WKM1 = TT*( ( D11*A( J, K-1 )-DCONJG( D12 )* + $ A( J, K ) ) / D ) + WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D ) * * Perform a rank-2 update of A(1:k-2,1:k-2) * DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*DCONJG( WK ) - - $ ( A( I, K-1 ) / D )*DCONJG( WKM1 ) + $ A( I, K )*DCONJG( WK ) - + $ A( I, K-1 )*DCONJG( WKM1 ) 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D - A( J, K-1 ) = WKM1 / D + A( J, K ) = WK + A( J, K-1 ) = WKM1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO ) * @@ -857,24 +857,24 @@ SUBROUTINE ZHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) ) - WKP1 = TT*( D22*A( J, K+1 )-DCONJG( D21 )* - $ A( J, K ) ) + WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D ) + WKP1 = TT*( ( D22*A( J, K+1 )-DCONJG( D21 )* + $ A( J, K ) ) / D ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N A( I, J ) = A( I, J ) - - $ ( A( I, K ) / D )*DCONJG( WK ) - - $ ( A( I, K+1 ) / D )*DCONJG( WKP1 ) + $ A( I, K )*DCONJG( WK ) - + $ A( I, K+1 )*DCONJG( WKP1 ) 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D - A( J, K+1 ) = WKP1 / D + A( J, K ) = WK + A( J, K+1 ) = WKP1 * (*) Make sure that diagonal element of pivot is real A( J, J ) = DCMPLX( DBLE( A( J, J ) ), ZERO ) * diff --git a/SRC/zhptrf.f b/SRC/zhptrf.f index 6558a1bb8..b7d22d2f0 100644 --- a/SRC/zhptrf.f +++ b/SRC/zhptrf.f @@ -389,13 +389,12 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) D11 = DBLE( AP( K+( K-1 )*K / 2 ) ) / D TT = ONE / ( D11*D22-ONE ) D12 = AP( K-1+( K-1 )*K / 2 ) / D - D = TT / D * DO 50 J = K - 2, 1, -1 - WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ DCONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) - WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12* - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = TT*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ DCONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) / D ) + WK = TT*( ( D22*AP( J+( K-1 )*K / 2 )-D12* + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*DCONJG( WK ) - @@ -603,14 +602,13 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) D22 = DBLE( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D TT = ONE / ( D11*D22-ONE ) D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D - D = TT / D * DO 100 J = K + 2, N - WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21* - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + WK = TT*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21* + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D ) + WKP1 = TT*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- $ DCONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / - $ 2 ) ) + $ 2 ) ) / D ) DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / diff --git a/SRC/zlahef.f b/SRC/zlahef.f index 039444a0d..3fdde8113 100644 --- a/SRC/zlahef.f +++ b/SRC/zlahef.f @@ -480,17 +480,17 @@ SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = * ( ( -1 ) ( D22 ) ) * -* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ), * ( ( -1 ) ( D22 ) ) * -* = ( conj(D21)*( D11 ) D21*( -1 ) ) -* ( ( -1 ) ( D22 ) ), -* * where D11 = d22/d21, * D22 = d11/conj(d21), -* D21 = T/d21, * T = 1/(D22*D11-1). * +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 or conj(d21) and then scaled by T. +* * (NOTE: No need to check for division by ZERO, * since that was ensured earlier in pivot search: * (a) d21 != 0, since in 2x2 pivot case(4) @@ -502,16 +502,16 @@ SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, D11 = W( K, KW ) / DCONJG( D21 ) D22 = W( K-1, KW-1 ) / D21 T = ONE / ( DBLE( D11*D22 )-ONE ) - D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = DCONJG( D21 )* - $ ( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( ( D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ DCONJG( D21 ) ) 20 CONTINUE END IF * @@ -827,17 +827,17 @@ SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = * ( ( -1 ) ( D22 ) ) * -* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ), * ( ( -1 ) ( D22 ) ) * -* = ( conj(D21)*( D11 ) D21*( -1 ) ) -* ( ( -1 ) ( D22 ) ), -* * where D11 = d22/d21, * D22 = d11/conj(d21), -* D21 = T/d21, * T = 1/(D22*D11-1). * +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 or conj(d21) and then scaled by T. +* * (NOTE: No need to check for division by ZERO, * since that was ensured earlier in pivot search: * (a) d21 != 0, since in 2x2 pivot case(4) @@ -849,16 +849,16 @@ SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / DCONJG( D21 ) T = ONE / ( DBLE( D11*D22 )-ONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = DCONJG( D21 )* - $ ( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ DCONJG( D21 ) ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/zlasyf.f b/SRC/zlasyf.f index 5363c390e..9411871fb 100644 --- a/SRC/zlasyf.f +++ b/SRC/zlasyf.f @@ -431,22 +431,24 @@ SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K-1, KW ) D11 = W( K, KW ) / D21 D22 = W( K-1, KW-1 ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 - A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) - A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) + A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) / + $ D21 ) + A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) / + $ D21 ) 20 CONTINUE END IF * @@ -710,22 +712,24 @@ SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * -* = D21 * ( ( D11 ) ( -1 ) ) -* ( ( -1 ) ( D22 ) ) +* T/d21 is not formed, since it overflows when d21 is +* subnormal: each entry of the product is divided by +* d21 and then scaled by T. * D21 = W( K+1, K ) D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N - A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) - A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) + A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) / + $ D21 ) + A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) / + $ D21 ) 80 CONTINUE END IF * diff --git a/SRC/zsptrf.f b/SRC/zsptrf.f index 4a5a7a50b..42a916cdb 100644 --- a/SRC/zsptrf.f +++ b/SRC/zsptrf.f @@ -375,13 +375,12 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12 D11 = AP( K+( K-1 )*K / 2 ) / D12 T = CONE / ( D11*D22-CONE ) - D12 = T / D12 * DO 50 J = K - 2, 1, -1 - WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )- - $ AP( J+( K-1 )*K / 2 ) ) - WK = D12*( D22*AP( J+( K-1 )*K / 2 )- - $ AP( J+( K-2 )*( K-1 ) / 2 ) ) + WKM1 = T*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )- + $ AP( J+( K-1 )*K / 2 ) ) / D12 ) + WK = T*( ( D22*AP( J+( K-1 )*K / 2 )- + $ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D12 ) DO 40 I = J, 1, -1 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) - $ AP( I+( K-1 )*K / 2 )*WK - @@ -576,13 +575,12 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21 D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * DO 100 J = K + 2, N - WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- - $ AP( J+K*( 2*N-K-1 ) / 2 ) ) - WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )- - $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) + WK = T*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )- + $ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D21 ) + WKP1 = T*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )- + $ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) / D21 ) DO 90 I = J, N AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )* $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) / diff --git a/SRC/zsytf2.f b/SRC/zsytf2.f index 63cb9709f..26ed82a08 100644 --- a/SRC/zsytf2.f +++ b/SRC/zsytf2.f @@ -394,11 +394,10 @@ SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D22 = A( K-1, K-1 ) / D12 D11 = A( K, K ) / D12 T = CONE / ( D11*D22-CONE ) - D12 = T / D12 * DO 30 J = K - 2, 1, -1 - WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) - WK = D12*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K-1 )*WKM1 @@ -569,11 +568,10 @@ SUBROUTINE ZSYTF2( UPLO, N, A, LDA, IPIV, INFO ) D11 = A( K+1, K+1 ) / D21 D22 = A( K, K ) / D21 T = CONE / ( D11*D22-CONE ) - D21 = T / D21 * DO 60 J = K + 2, N - WK = D21*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) DO 50 I = J, N A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K+1 )*WKP1 diff --git a/SRC/zsytf2_rk.f b/SRC/zsytf2_rk.f index 9d549adc5..66bac797a 100644 --- a/SRC/zsytf2_rk.f +++ b/SRC/zsytf2_rk.f @@ -582,18 +582,18 @@ SUBROUTINE ZSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -898,22 +898,22 @@ SUBROUTINE ZSYTF2_RK( UPLO, N, A, LDA, E, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/SRC/zsytf2_rook.f b/SRC/zsytf2_rook.f index 9038e0bbd..5ca176196 100644 --- a/SRC/zsytf2_rook.f +++ b/SRC/zsytf2_rook.f @@ -502,18 +502,18 @@ SUBROUTINE ZSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 30 J = K - 2, 1, -1 * - WKM1 = T*( D11*A( J, K-1 )-A( J, K ) ) - WK = T*( D22*A( J, K )-A( J, K-1 ) ) + WKM1 = T*( ( D11*A( J, K-1 )-A( J, K ) ) / D12 ) + WK = T*( ( D22*A( J, K )-A( J, K-1 ) ) / D12 ) * DO 20 I = J, 1, -1 - A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK - - $ ( A( I, K-1 ) / D12 )*WKM1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K-1 )*WKM1 20 CONTINUE * * Store U(k) and U(k-1) in cols k and k-1 for row J * - A( J, K ) = WK / D12 - A( J, K-1 ) = WKM1 / D12 + A( J, K ) = WK + A( J, K-1 ) = WKM1 * 30 CONTINUE * @@ -776,22 +776,22 @@ SUBROUTINE ZSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO ) * DO 60 J = K + 2, N * -* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J +* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J * - WK = T*( D11*A( J, K )-A( J, K+1 ) ) - WKP1 = T*( D22*A( J, K+1 )-A( J, K ) ) + WK = T*( ( D11*A( J, K )-A( J, K+1 ) ) / D21 ) + WKP1 = T*( ( D22*A( J, K+1 )-A( J, K ) ) / D21 ) * * Perform a rank-2 update of A(k+2:n,k+2:n) * DO 50 I = J, N - A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK - - $ ( A( I, K+1 ) / D21 )*WKP1 + A( I, J ) = A( I, J ) - A( I, K )*WK - + $ A( I, K+1 )*WKP1 50 CONTINUE * * Store L(k) and L(k+1) in cols k and k+1 for row J * - A( J, K ) = WK / D21 - A( J, K+1 ) = WKP1 / D21 + A( J, K ) = WK + A( J, K+1 ) = WKP1 * 60 CONTINUE * diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index b04a3f796..880effbeb 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -300,7 +300,7 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) END IF WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' ) IF( SORD ) THEN - WRITE( IOUNIT, FMT = 9972 ) + WRITE( IOUNIT, FMT = 7972 ) ELSE WRITE( IOUNIT, FMT = 9971 ) END IF @@ -913,6 +913,21 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) $ 'TRF, no test ratios are computed)' ) * * SSY, SSR, SSP, CHE, CHR, CHP matrix types +* +* +* SSY matrix types +* + 7972 FORMAT( 4X, '1. Diagonal', 24X, + $ '6. Last n/2 rows and columns zero', / 4X, + $ '2. Random, CNDNUM = 2', 14X, + $ '7. Random, CNDNUM = sqrt(0.1/EPS)', / 4X, + $ '3. First row and column zero', 7X, + $ '8. Random, CNDNUM = 0.1/EPS', / 4X, + $ '4. Last row and column zero', 8X, + $ '9. Scaled near underflow', / 4X, + $ '5. Middle row and column zero', 5X, + $ '10. Scaled near overflow', / 39X, + $ '11. Subnormal 2 by 2 pivot block' ) * 9972 FORMAT( 4X, '1. Diagonal', 24X, $ '6. Last n/2 rows and columns zero', / 4X, diff --git a/TESTING/LIN/dchkaa.F b/TESTING/LIN/dchkaa.F index d0f7dbcb4..ed5c7ab90 100644 --- a/TESTING/LIN/dchkaa.F +++ b/TESTING/LIN/dchkaa.F @@ -671,7 +671,7 @@ PROGRAM DCHKAA * SY: symmetric indefinite matrices, * with partial (Bunch-Kaufman) pivoting algorithm * - NTYPES = 10 + NTYPES = 11 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/dchksy.f b/TESTING/LIN/dchksy.f index 4d9789e74..cf4fc5e4d 100644 --- a/TESTING/LIN/dchksy.f +++ b/TESTING/LIN/dchksy.f @@ -177,6 +177,7 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, LOGICAL TSTERR INTEGER NMAX, NN, NNB, NNS, NOUT DOUBLE PRECISION THRESH + DOUBLE PRECISION SUBNRM * .. * .. Array Arguments .. LOGICAL DOTYPE( * ) @@ -190,8 +191,10 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) + DOUBLE PRECISION FOUR + PARAMETER ( FOUR = 4.0D+0 ) INTEGER NTYPES - PARAMETER ( NTYPES = 10 ) + PARAMETER ( NTYPES = 11 ) INTEGER NTESTS PARAMETER ( NTESTS = 9 ) * .. @@ -210,13 +213,16 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH DOUBLE PRECISION DGET06, DLANSY EXTERNAL DGET06, DLANSY + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY, $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05, - $ DSYCON, DSYRFS, DSYT01, DSYTRF, + $ DSCAL, DSYCON, DSYRFS, DSYT01, DSYTRF, $ DSYTRI2, DSYTRS, DSYTRS2, XLAENV * .. * .. Intrinsic Functions .. @@ -387,6 +393,32 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, IZERO = 0 END IF * +* Type 11: scale two adjacent rows and columns into the +* subnormal range and give them a 2 by 2 pivot block whose +* off-diagonal entry is four times its diagonal, at the +* end the factorization starts from. Inverting that pivot +* through the reciprocal of the off-diagonal entry +* overflows. +* + IF( IMAT.EQ.11 .AND. N.GE.2 ) THEN + SUBNRM = DLAMCH( 'Safe minimum' ) / 512 + IF( IUPLO.EQ.1 ) THEN + I1 = N - 1 + ELSE + I1 = 1 + END IF + I2 = I1 + 1 + CALL DSCAL( N, SUBNRM, A( I1 ), LDA ) + CALL DSCAL( N, SUBNRM, A( I2 ), LDA ) + CALL DSCAL( N, SUBNRM, A( ( I1-1 )*LDA+1 ), 1 ) + CALL DSCAL( N, SUBNRM, A( ( I2-1 )*LDA+1 ), 1 ) + A( ( I1-1 )*LDA+I1 ) = SUBNRM + A( ( I2-1 )*LDA+I2 ) = SUBNRM + A( ( I2-1 )*LDA+I1 ) = FOUR*SUBNRM + A( ( I1-1 )*LDA+I2 ) = FOUR*SUBNRM + IZERO = 0 + END IF +* * End generate the test matrix A. * * Do for each value of NB in NBVAL @@ -440,7 +472,7 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * * Set the condition estimate flag if the INFO is not 0. * - IF( INFO.NE.0 ) THEN + IF( INFO.NE.0 .OR. IMAT.EQ.11 ) THEN TRFCON = .TRUE. ELSE TRFCON = .FALSE. @@ -485,7 +517,8 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * the threshold. * DO 110 K = 1, NT - IF( RESULT( K ).GE.THRESH ) THEN + IF( RESULT( K ).GE.THRESH .OR. + $ DISNAN( RESULT( K ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, @@ -624,6 +657,8 @@ SUBROUTINE DCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE + IF( IMAT.EQ.11 ) + $ GO TO 150 ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK ) SRNAMT = 'DSYCON' CALL DSYCON( UPLO, N, AFAC, LDA, IWORK, ANORM, RCOND, diff --git a/TESTING/LIN/schkaa.F b/TESTING/LIN/schkaa.F index dd3f8de4b..bb5b0af2c 100644 --- a/TESTING/LIN/schkaa.F +++ b/TESTING/LIN/schkaa.F @@ -667,7 +667,7 @@ PROGRAM SCHKAA * SY: symmetric indefinite matrices, * with partial (Bunch-Kaufman) pivoting algorithm * - NTYPES = 10 + NTYPES = 11 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/schksy.f b/TESTING/LIN/schksy.f index a8de72ca6..12f416df9 100644 --- a/TESTING/LIN/schksy.f +++ b/TESTING/LIN/schksy.f @@ -177,6 +177,7 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, LOGICAL TSTERR INTEGER NMAX, NN, NNB, NNS, NOUT REAL THRESH + REAL SUBNRM * .. * .. Array Arguments .. LOGICAL DOTYPE( * ) @@ -190,8 +191,10 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) + REAL FOUR + PARAMETER ( FOUR = 4.0E+0 ) INTEGER NTYPES - PARAMETER ( NTYPES = 10 ) + PARAMETER ( NTYPES = 11 ) INTEGER NTESTS PARAMETER ( NTESTS = 9 ) * .. @@ -210,14 +213,17 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH REAL SGET06, SLANSY EXTERNAL SGET06, SLANSY + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRSY, SGET04, SLACPY, $ SLARHS, SLATB4, SLATMS, SPOT02, SPOT03, SPOT05, - $ SSYCON, SSYRFS, SSYT01, SSYTRF, SSYTRI2, - $ SSYTRS, SSYTRS2, XLAENV + $ SSCAL, SSYCON, SSYRFS, SSYT01, SSYTRF, + $ SSYTRI2, SSYTRS, SSYTRS2, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -386,6 +392,32 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, IZERO = 0 END IF * +* Type 11: scale two adjacent rows and columns into the +* subnormal range and give them a 2 by 2 pivot block whose +* off-diagonal entry is four times its diagonal, at the +* end the factorization starts from. Inverting that pivot +* through the reciprocal of the off-diagonal entry +* overflows. +* + IF( IMAT.EQ.11 .AND. N.GE.2 ) THEN + SUBNRM = SLAMCH( 'Safe minimum' ) / 512 + IF( IUPLO.EQ.1 ) THEN + I1 = N - 1 + ELSE + I1 = 1 + END IF + I2 = I1 + 1 + CALL SSCAL( N, SUBNRM, A( I1 ), LDA ) + CALL SSCAL( N, SUBNRM, A( I2 ), LDA ) + CALL SSCAL( N, SUBNRM, A( ( I1-1 )*LDA+1 ), 1 ) + CALL SSCAL( N, SUBNRM, A( ( I2-1 )*LDA+1 ), 1 ) + A( ( I1-1 )*LDA+I1 ) = SUBNRM + A( ( I2-1 )*LDA+I2 ) = SUBNRM + A( ( I2-1 )*LDA+I1 ) = FOUR*SUBNRM + A( ( I1-1 )*LDA+I2 ) = FOUR*SUBNRM + IZERO = 0 + END IF +* * End generate the test matrix A. * * @@ -440,7 +472,7 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * * Set the condition estimate flag if the INFO is not 0. * - IF( INFO.NE.0 ) THEN + IF( INFO.NE.0 .OR. IMAT.EQ.11 ) THEN TRFCON = .TRUE. ELSE TRFCON = .FALSE. @@ -485,7 +517,8 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * the threshold. * DO 110 K = 1, NT - IF( RESULT( K ).GE.THRESH ) THEN + IF( RESULT( K ).GE.THRESH .OR. + $ SISNAN( RESULT( K ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K, @@ -623,6 +656,8 @@ SUBROUTINE SCHKSY( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL, * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE + IF( IMAT.EQ.11 ) + $ GO TO 150 ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK ) SRNAMT = 'SSYCON' CALL SSYCON( UPLO, N, AFAC, LDA, IWORK, ANORM, RCOND, diff --git a/TESTING/dtest.in b/TESTING/dtest.in index cde62db50..34a6dd14b 100644 --- a/TESTING/dtest.in +++ b/TESTING/dtest.in @@ -22,7 +22,7 @@ DPS 9 List types on next line if 0 < NTYPES < 9 DPP 9 List types on next line if 0 < NTYPES < 9 DPB 8 List types on next line if 0 < NTYPES < 8 DPT 12 List types on next line if 0 < NTYPES < 12 -DSY 10 List types on next line if 0 < NTYPES < 10 +DSY 11 List types on next line if 0 < NTYPES < 11 DSR 10 List types on next line if 0 < NTYPES < 10 DSK 10 List types on next line if 0 < NTYPES < 10 DSA 10 List types on next line if 0 < NTYPES < 10 diff --git a/TESTING/stest.in b/TESTING/stest.in index abfd639fd..b95b4ed32 100644 --- a/TESTING/stest.in +++ b/TESTING/stest.in @@ -22,7 +22,7 @@ SPS 9 List types on next line if 0 < NTYPES < 9 SPP 9 List types on next line if 0 < NTYPES < 9 SPB 8 List types on next line if 0 < NTYPES < 8 SPT 12 List types on next line if 0 < NTYPES < 12 -SSY 10 List types on next line if 0 < NTYPES < 10 +SSY 11 List types on next line if 0 < NTYPES < 11 SSR 10 List types on next line if 0 < NTYPES < 10 SSK 10 List types on next line if 0 < NTYPES < 10 SSA 10 List types on next line if 0 < NTYPES < 10