Update LAPACK to 3.8.0

This commit is contained in:
martin
2017-11-23 18:13:35 +01:00
parent b18730f9e1
commit 3be5c3d343
2092 changed files with 45721 additions and 23035 deletions
+59 -72
View File
@@ -19,11 +19,11 @@
* ===========
*
* SUBROUTINE CLASYF_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 REAL 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 complexSYcomputational
*
* =====================================================================
SUBROUTINE CLASYF_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.0E+0, ONE = 1.0E+0 )
*
* .. Local Scalars ..
INTEGER J, K, K1, I1, I2
INTEGER J, K, K1, I1, I2, MJ
COMPLEX PIV, ALPHA
* ..
* .. External Functions ..
@@ -185,14 +175,14 @@
EXTERNAL LSAME, ILAENV, ICAMAX
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
EXTERNAL CAXPY, CGEMV, CSCAL, CCOPY, CSWAP, CLASET,
$ 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
*
* 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)
* Only need to compute T(J, J)
*
MJ = 1
ELSE
MJ = M-J+1
END IF
*
* 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 CGEMV( 'No transpose', M-J+1, J-K1,
CALL CGEMV( '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 CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
CALL CCOPY( 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 CAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
CALL CAXPY( 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 = ICAMAX( 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 CSWAP( 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 CSWAP( 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 CCOPY( 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 CLASET( '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
*
* 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)
* Only need to compute T(J, J)
*
MJ = 1
ELSE
MJ = M-J+1
END IF
*
* 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 CGEMV( 'No transpose', M-J+1, J-K1,
CALL CGEMV( '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 CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
CALL CCOPY( 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 CAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
CALL CAXPY( 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 = ICAMAX( 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 CSWAP( 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 CSWAP( 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 CCOPY( 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 CLASET( '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