version 1.1, 2017/06/17 11:02:50
|
version 1.4, 2018/05/29 07:18:01
|
Line 19
|
Line 19
|
* =========== |
* =========== |
* |
* |
* SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, |
* SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, |
* H, LDH, WORK, INFO ) |
* H, LDH, WORK ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER UPLO |
* CHARACTER UPLO |
* INTEGER J1, M, NB, LDA, LDH, INFO |
* INTEGER J1, M, NB, LDA, LDH |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
* INTEGER IPIV( * ) |
* INTEGER IPIV( * ) |
Line 99
|
Line 99
|
*> \param[in] LDA |
*> \param[in] LDA |
*> \verbatim |
*> \verbatim |
*> LDA is INTEGER |
*> 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 |
*> \endverbatim |
*> |
*> |
*> \param[out] IPIV |
*> \param[out] IPIV |
*> \verbatim |
*> \verbatim |
*> IPIV is INTEGER array, dimension (N) |
*> IPIV is INTEGER array, dimension (M) |
*> Details of the row and column interchanges, |
*> Details of the row and column interchanges, |
*> the row and column k were interchanged with the row and |
*> the row and column k were interchanged with the row and |
*> column IPIV(k). |
*> column IPIV(k). |
Line 127
|
Line 127
|
*> WORK is DOUBLE PRECISION workspace, dimension (M). |
*> WORK is DOUBLE PRECISION workspace, dimension (M). |
*> \endverbatim |
*> \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: |
* Authors: |
* ======== |
* ======== |
Line 146
|
Line 136
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date December 2016 |
*> \date November 2017 |
* |
* |
*> \ingroup doubleSYcomputational |
*> \ingroup doubleSYcomputational |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, |
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, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* December 2016 |
* November 2017 |
* |
* |
IMPLICIT NONE |
IMPLICIT NONE |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER UPLO |
CHARACTER UPLO |
INTEGER M, NB, J1, LDA, LDH, INFO |
INTEGER M, NB, J1, LDA, LDH |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
INTEGER IPIV( * ) |
INTEGER IPIV( * ) |
Line 176
|
Line 166
|
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
* |
* |
* .. Local Scalars .. |
* .. Local Scalars .. |
INTEGER J, K, K1, I1, I2 |
INTEGER J, K, K1, I1, I2, MJ |
DOUBLE PRECISION PIV, ALPHA |
DOUBLE PRECISION PIV, ALPHA |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
Line 185
|
Line 175
|
EXTERNAL LSAME, ILAENV, IDAMAX |
EXTERNAL LSAME, ILAENV, IDAMAX |
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL XERBLA |
EXTERNAL DGEMV, DAXPY, DCOPY, DSWAP, DSCAL, DLASET, |
|
$ XERBLA |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC MAX |
INTRINSIC MAX |
* .. |
* .. |
* .. Executable Statements .. |
* .. Executable Statements .. |
* |
* |
INFO = 0 |
|
J = 1 |
J = 1 |
* |
* |
* K1 is the first column of the panel to be factorized |
* K1 is the first column of the panel to be factorized |
Line 216
|
Line 206
|
* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, |
* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, |
* |
* |
K = J1+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), |
* H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J), |
* where H(J:N, J) has been initialized to be A(J, J:N) |
* where H(J:M, J) has been initialized to be A(J, J:M) |
* |
* |
IF( K.GT.2 ) THEN |
IF( K.GT.2 ) THEN |
* |
* |
Line 228
|
Line 226
|
* > for the rest of the columns, K is J+1, skipping only the |
* > for the rest of the columns, K is J+1, skipping only the |
* first column |
* first column |
* |
* |
CALL DGEMV( 'No transpose', M-J+1, J-K1, |
CALL DGEMV( 'No transpose', MJ, J-K1, |
$ -ONE, H( J, K1 ), LDH, |
$ -ONE, H( J, K1 ), LDH, |
$ A( 1, J ), 1, |
$ A( 1, J ), 1, |
$ ONE, H( J, J ), 1 ) |
$ ONE, H( J, J ), 1 ) |
END IF |
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 |
IF( J.GT.K1 ) THEN |
* |
* |
* Compute WORK := WORK - L(J-1, J:N) * T(J-1,J), |
* 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:N) stores U(J-1, J:N) |
* 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 ) |
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 |
END IF |
* |
* |
* Set A(J, J) = T(J, J) |
* Set A(J, J) = T(J, J) |
Line 253
|
Line 251
|
* |
* |
IF( J.LT.M ) THEN |
IF( J.LT.M ) THEN |
* |
* |
* Compute WORK(2:N) = T(J, J) L(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):N) stores U(J, (J+1):N) |
* 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 |
IF( K.GT.1 ) THEN |
ALPHA = -A( K, J ) |
ALPHA = -A( K, J ) |
Line 262
|
Line 260
|
$ WORK( 2 ), 1 ) |
$ WORK( 2 ), 1 ) |
ENDIF |
ENDIF |
* |
* |
* Find max(|WORK(2:n)|) |
* Find max(|WORK(2:M)|) |
* |
* |
I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 |
I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 |
PIV = WORK( I2 ) |
PIV = WORK( I2 ) |
Line 277
|
Line 275
|
WORK( I2 ) = WORK( I1 ) |
WORK( I2 ) = WORK( I1 ) |
WORK( I1 ) = PIV |
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 |
I1 = I1+J-1 |
I2 = I2+J-1 |
I2 = I2+J-1 |
CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, |
CALL DSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, |
$ A( J1+I1, I2 ), 1 ) |
$ 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, |
CALL DSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, |
$ A( J1+I2-1, I2+1 ), LDA ) |
$ A( J1+I2-1, I2+1 ), LDA ) |
Line 315
|
Line 313
|
* Set A(J, J+1) = T(J, J+1) |
* Set A(J, J+1) = T(J, J+1) |
* |
* |
A( K, J+1 ) = WORK( 2 ) |
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 |
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, |
CALL DCOPY( M-J, A( K+1, J+1 ), LDA, |
$ H( J+1, J+1 ), 1 ) |
$ H( J+1, J+1 ), 1 ) |
END IF |
END IF |
* |
* |
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, 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:N, J) = L(J+2:N, 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 |
IF( A( K, J+1 ).NE.ZERO ) THEN |
ALPHA = ONE / A( K, J+1 ) |
ALPHA = ONE / A( K, J+1 ) |
Line 341
|
Line 333
|
CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO, |
CALL DLASET( 'Full', 1, M-J-1, ZERO, ZERO, |
$ A( K, J+2 ), LDA) |
$ A( K, J+2 ), LDA) |
END IF |
END IF |
ELSE |
|
IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN |
|
INFO = J |
|
END IF |
|
END IF |
END IF |
J = J + 1 |
J = J + 1 |
GO TO 10 |
GO TO 10 |
Line 366
|
Line 354
|
* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, |
* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, |
* |
* |
K = J1+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, |
* H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T, |
* where H(J:N, J) has been initialized to be A(J:N, J) |
* where H(J:M, J) has been initialized to be A(J:M, J) |
* |
* |
IF( K.GT.2 ) THEN |
IF( K.GT.2 ) THEN |
* |
* |
Line 378
|
Line 374
|
* > for the rest of the columns, K is J+1, skipping only the |
* > for the rest of the columns, K is J+1, skipping only the |
* first column |
* first column |
* |
* |
CALL DGEMV( 'No transpose', M-J+1, J-K1, |
CALL DGEMV( 'No transpose', MJ, J-K1, |
$ -ONE, H( J, K1 ), LDH, |
$ -ONE, H( J, K1 ), LDH, |
$ A( J, 1 ), LDA, |
$ A( J, 1 ), LDA, |
$ ONE, H( J, J ), 1 ) |
$ ONE, H( J, J ), 1 ) |
END IF |
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 |
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) |
* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1) |
* |
* |
ALPHA = -A( J, K-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 |
END IF |
* |
* |
* Set A(J, J) = T(J, J) |
* Set A(J, J) = T(J, J) |
Line 403
|
Line 399
|
* |
* |
IF( J.LT.M ) THEN |
IF( J.LT.M ) THEN |
* |
* |
* Compute WORK(2:N) = T(J, J) 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):N, J-1) = L((J+1):N, 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 |
IF( K.GT.1 ) THEN |
ALPHA = -A( J, K ) |
ALPHA = -A( J, K ) |
Line 412
|
Line 408
|
$ WORK( 2 ), 1 ) |
$ WORK( 2 ), 1 ) |
ENDIF |
ENDIF |
* |
* |
* Find max(|WORK(2:n)|) |
* Find max(|WORK(2:M)|) |
* |
* |
I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 |
I2 = IDAMAX( M-J, WORK( 2 ), 1 ) + 1 |
PIV = WORK( I2 ) |
PIV = WORK( I2 ) |
Line 427
|
Line 423
|
WORK( I2 ) = WORK( I1 ) |
WORK( I2 ) = WORK( I1 ) |
WORK( I1 ) = PIV |
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 |
I1 = I1+J-1 |
I2 = I2+J-1 |
I2 = I2+J-1 |
CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, |
CALL DSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, |
$ A( I2, J1+I1 ), LDA ) |
$ 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, |
CALL DSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, |
$ A( I2+1, J1+I2-1 ), 1 ) |
$ A( I2+1, J1+I2-1 ), 1 ) |
Line 465
|
Line 461
|
* Set A(J+1, J) = T(J+1, J) |
* Set A(J+1, J) = T(J+1, J) |
* |
* |
A( J+1, K ) = WORK( 2 ) |
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 |
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, |
CALL DCOPY( M-J, A( J+1, K+1 ), 1, |
$ H( J+1, J+1 ), 1 ) |
$ H( J+1, J+1 ), 1 ) |
END IF |
END IF |
* |
* |
* Compute L(J+2, J+1) = WORK( 3:N ) / T(J, 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:N, J) = L(J+2:N, 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 |
IF( A( J+1, K ).NE.ZERO ) THEN |
ALPHA = ONE / A( J+1, K ) |
ALPHA = ONE / A( J+1, K ) |
Line 490
|
Line 481
|
CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO, |
CALL DLASET( 'Full', M-J-1, 1, ZERO, ZERO, |
$ A( J+2, K ), LDA ) |
$ A( J+2, K ), LDA ) |
END IF |
END IF |
ELSE |
|
IF( (A( J, K ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN |
|
INFO = J |
|
END IF |
|
END IF |
END IF |
J = J + 1 |
J = J + 1 |
GO TO 30 |
GO TO 30 |