--- rpl/lapack/lapack/dlasyf_aa.f 2017/06/17 11:06:27 1.2 +++ rpl/lapack/lapack/dlasyf_aa.f 2018/05/29 06:55:20 1.3 @@ -19,11 +19,11 @@ * =========== * * SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, -* H, LDH, WORK, INFO ) +* H, LDH, WORK ) * * .. Scalar Arguments .. * CHARACTER UPLO -* INTEGER J1, M, NB, LDA, LDH, INFO +* INTEGER J1, M, NB, LDA, LDH * .. * .. Array Arguments .. * INTEGER IPIV( * ) @@ -99,12 +99,12 @@ *> \param[in] LDA *> \verbatim *> LDA is INTEGER -*> The leading dimension of the array A. LDA >= max(1,N). +*> The leading dimension of the array A. LDA >= max(1,M). *> \endverbatim *> *> \param[out] IPIV *> \verbatim -*> IPIV is INTEGER array, dimension (N) +*> IPIV is INTEGER array, dimension (M) *> Details of the row and column interchanges, *> the row and column k were interchanged with the row and *> column IPIV(k). @@ -127,16 +127,6 @@ *> WORK is DOUBLE PRECISION workspace, dimension (M). *> \endverbatim *> -*> \param[out] INFO -*> \verbatim -*> INFO is INTEGER -*> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> > 0: if INFO = i, D(i,i) is exactly zero. The factorization -*> has been completed, but the block diagonal matrix D is -*> exactly singular, and division by zero will occur if it -*> is used to solve a system of equations. -*> \endverbatim * * Authors: * ======== @@ -146,24 +136,24 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2017 * *> \ingroup doubleSYcomputational * * ===================================================================== SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, - $ H, LDH, WORK, INFO ) + $ H, LDH, WORK ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.8.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2017 * IMPLICIT NONE * * .. Scalar Arguments .. CHARACTER UPLO - INTEGER M, NB, J1, LDA, LDH, INFO + INTEGER M, NB, J1, LDA, LDH * .. * .. Array Arguments .. INTEGER IPIV( * ) @@ -176,7 +166,7 @@ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * * .. Local Scalars .. - INTEGER J, K, K1, I1, I2 + INTEGER J, K, K1, I1, I2, MJ DOUBLE PRECISION PIV, ALPHA * .. * .. External Functions .. @@ -185,14 +175,14 @@ EXTERNAL LSAME, ILAENV, IDAMAX * .. * .. External Subroutines .. - EXTERNAL XERBLA + EXTERNAL DGEMV, DAXPY, DCOPY, DSWAP, DSCAL, DLASET, + $ XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * - INFO = 0 J = 1 * * K1 is the first column of the panel to be factorized @@ -216,9 +206,17 @@ * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, * K = J1+J-1 + IF( J.EQ.M ) THEN +* +* Only need to compute T(J, J) +* + MJ = 1 + ELSE + MJ = M-J+1 + END IF * -* H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J), -* where H(J:N, J) has been initialized to be A(J, J:N) +* H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J), +* where H(J:M, J) has been initialized to be A(J, J:M) * IF( K.GT.2 ) THEN * @@ -228,23 +226,23 @@ * > for the rest of the columns, K is J+1, skipping only the * first column * - CALL DGEMV( 'No transpose', M-J+1, J-K1, + CALL DGEMV( 'No transpose', MJ, J-K1, $ -ONE, H( J, K1 ), LDH, $ A( 1, J ), 1, $ ONE, H( J, J ), 1 ) END IF * -* Copy H(i:n, i) into WORK +* Copy H(i:M, i) into WORK * - CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) + CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 ) * IF( J.GT.K1 ) THEN * -* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J), -* where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N) +* Compute WORK := WORK - L(J-1, J:M) * T(J-1,J), +* where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M) * ALPHA = -A( K-1, J ) - CALL DAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 ) + CALL DAXPY( MJ, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 ) END IF * * Set A(J, J) = T(J, J) @@ -253,8 +251,8 @@ * IF( J.LT.M ) THEN * -* Compute WORK(2:N) = T(J, J) L(J, (J+1):N) -* where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N) +* Compute WORK(2:M) = T(J, J) L(J, (J+1):M) +* where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M) * IF( K.GT.1 ) THEN ALPHA = -A( K, J ) @@ -262,7 +260,7 @@ $ WORK( 2 ), 1 ) ENDIF * -* Find max(|WORK(2:n)|) +* Find max(|WORK(2:M)|) * I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 PIV = WORK( I2 ) @@ -277,14 +275,14 @@ WORK( I2 ) = WORK( I1 ) WORK( I1 ) = PIV * -* Swap A(I1, I1+1:N) with A(I1+1:N, I2) +* Swap A(I1, I1+1:M) with A(I1+1:M, I2) * I1 = I1+J-1 I2 = I2+J-1 CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, $ A( J1+I1, I2 ), 1 ) * -* Swap A(I1, I2+1:N) with A(I2, I2+1:N) +* Swap A(I1, I2+1:M) with A(I2, I2+1:M) * CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, $ A( J1+I2-1, I2+1 ), LDA ) @@ -315,23 +313,17 @@ * Set A(J, J+1) = T(J, J+1) * A( K, J+1 ) = WORK( 2 ) - IF( (A( K, J ).EQ.ZERO ) .AND. - $ ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN - IF(INFO .EQ. 0) THEN - INFO = J - ENDIF - END IF * IF( J.LT.NB ) THEN * -* Copy A(J+1:N, J+1) into H(J:N, J), +* Copy A(J+1:M, J+1) into H(J:M, J), * CALL DCOPY( M-J, A( K+1, J+1 ), LDA, $ H( J+1, J+1 ), 1 ) END IF * -* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1), -* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1) +* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), +* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * IF( A( K, J+1 ).NE.ZERO ) THEN ALPHA = ONE / A( K, J+1 ) @@ -341,10 +333,6 @@ CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO, $ A( K, J+2 ), LDA) END IF - ELSE - IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN - INFO = J - END IF END IF J = J + 1 GO TO 10 @@ -366,9 +354,17 @@ * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, * K = J1+J-1 + IF( J.EQ.M ) THEN +* +* Only need to compute T(J, J) +* + MJ = 1 + ELSE + MJ = M-J+1 + END IF * -* H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T, -* where H(J:N, J) has been initialized to be A(J:N, J) +* H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T, +* where H(J:M, J) has been initialized to be A(J:M, J) * IF( K.GT.2 ) THEN * @@ -378,23 +374,23 @@ * > for the rest of the columns, K is J+1, skipping only the * first column * - CALL DGEMV( 'No transpose', M-J+1, J-K1, + CALL DGEMV( 'No transpose', MJ, J-K1, $ -ONE, H( J, K1 ), LDH, $ A( J, 1 ), LDA, $ ONE, H( J, J ), 1 ) END IF * -* Copy H(J:N, J) into WORK +* Copy H(J:M, J) into WORK * - CALL DCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 ) + CALL DCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 ) * IF( J.GT.K1 ) THEN * -* Compute WORK := WORK - L(J:N, J-1) * T(J-1,J), +* Compute WORK := WORK - L(J:M, J-1) * T(J-1,J), * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1) * ALPHA = -A( J, K-1 ) - CALL DAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 ) + CALL DAXPY( MJ, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 ) END IF * * Set A(J, J) = T(J, J) @@ -403,8 +399,8 @@ * IF( J.LT.M ) THEN * -* Compute WORK(2:N) = T(J, J) L((J+1):N, J) -* where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J) +* Compute WORK(2:M) = T(J, J) L((J+1):M, J) +* where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J) * IF( K.GT.1 ) THEN ALPHA = -A( J, K ) @@ -412,7 +408,7 @@ $ WORK( 2 ), 1 ) ENDIF * -* Find max(|WORK(2:n)|) +* Find max(|WORK(2:M)|) * I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 PIV = WORK( I2 ) @@ -427,14 +423,14 @@ WORK( I2 ) = WORK( I1 ) WORK( I1 ) = PIV * -* Swap A(I1+1:N, I1) with A(I2, I1+1:N) +* Swap A(I1+1:M, I1) with A(I2, I1+1:M) * I1 = I1+J-1 I2 = I2+J-1 CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, $ A( I2, J1+I1 ), LDA ) * -* Swap A(I2+1:N, I1) with A(I2+1:N, I2) +* Swap A(I2+1:M, I1) with A(I2+1:M, I2) * CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, $ A( I2+1, J1+I2-1 ), 1 ) @@ -465,22 +461,17 @@ * Set A(J+1, J) = T(J+1, J) * A( J+1, K ) = WORK( 2 ) - IF( (A( J, K ).EQ.ZERO) .AND. - $ ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN - IF (INFO .EQ. 0) - $ INFO = J - END IF * IF( J.LT.NB ) THEN * -* Copy A(J+1:N, J+1) into H(J+1:N, J), +* Copy A(J+1:M, J+1) into H(J+1:M, J), * CALL DCOPY( M-J, A( J+1, K+1 ), 1, $ H( J+1, J+1 ), 1 ) END IF * -* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1), -* where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1) +* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), +* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) * IF( A( J+1, K ).NE.ZERO ) THEN ALPHA = ONE / A( J+1, K ) @@ -490,10 +481,6 @@ CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO, $ A( J+2, K ), LDA ) END IF - ELSE - IF( (A( J, K ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN - INFO = J - END IF END IF J = J + 1 GO TO 30