--- rpl/lapack/lapack/dlasyf.f 2012/08/22 09:48:21 1.11
+++ rpl/lapack/lapack/dlasyf.f 2018/05/29 07:18:01 1.18
@@ -1,25 +1,25 @@
-*> \brief \b DLASYF
+*> \brief \b DLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal pivoting method.
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DLASYF + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download DLASYF + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
-*
+*
* .. Scalar Arguments ..
* CHARACTER UPLO
* INTEGER INFO, KB, LDA, LDW, N, NB
@@ -28,7 +28,7 @@
* INTEGER IPIV( * )
* DOUBLE PRECISION A( LDA, * ), W( LDW, * )
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -109,16 +109,26 @@
*> \verbatim
*> IPIV is INTEGER array, dimension (N)
*> Details of the interchanges and the block structure of D.
-*> If UPLO = 'U', only the last KB elements of IPIV are set;
-*> if UPLO = 'L', only the first KB elements are set.
*>
-*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
-*> interchanged and D(k,k) is a 1-by-1 diagonal block.
-*> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
-*> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
-*> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
-*> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
-*> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
+*> If UPLO = 'U':
+*> Only the last KB elements of IPIV are set.
+*>
+*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
+*> interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
+*> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
+*> is a 2-by-2 diagonal block.
+*>
+*> If UPLO = 'L':
+*> Only the first KB elements of IPIV are set.
+*>
+*> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
+*> interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
+*> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
+*> is a 2-by-2 diagonal block.
*> \endverbatim
*>
*> \param[out] W
@@ -144,22 +154,32 @@
* Authors:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
-*> \date November 2011
+*> \date November 2013
*
*> \ingroup doubleSYcomputational
*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2013, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*> \endverbatim
+*
* =====================================================================
SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
*
-* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK computational routine (version 3.5.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2011
+* November 2013
*
* .. Scalar Arguments ..
CHARACTER UPLO
@@ -237,7 +257,8 @@
ABSAKK = ABS( W( K, KW ) )
*
* IMAX is the row-index of the largest off-diagonal element in
-* column K, and COLMAX is its absolute value
+* column K, and COLMAX is its absolute value.
+* Determine both COLMAX and IMAX.
*
IF( K.GT.1 ) THEN
IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
@@ -248,7 +269,7 @@
*
IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
*
-* Column K is zero: set INFO and continue
+* Column K is zero or underflow: set INFO and continue
*
IF( INFO.EQ.0 )
$ INFO = K
@@ -293,7 +314,7 @@
*
KP = IMAX
*
-* copy column KW-1 of W to column KW
+* copy column KW-1 of W to column KW of W
*
CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
ELSE
@@ -306,59 +327,117 @@
END IF
END IF
*
+* ============================================================
+*
+* KK is the column of A where pivoting step stopped
+*
KK = K - KSTEP + 1
+*
+* KKW is the column of W which corresponds to column KK of A
+*
KKW = NB + KK - N
*
-* Updated column KP is already stored in column KKW of W
+* Interchange rows and columns KP and KK.
+* Updated column KP is already stored in column KKW of W.
*
IF( KP.NE.KK ) THEN
*
-* Copy non-updated column KK to column KP
+* Copy non-updated column KK to column KP of submatrix A
+* at step K. No need to copy element into column K
+* (or K and K-1 for 2-by-2 pivot) of A, since these columns
+* will be later overwritten.
*
- A( KP, K ) = A( KK, K )
- CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
+ A( KP, KP ) = A( KK, KK )
+ CALL DCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
$ LDA )
- CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
+ IF( KP.GT.1 )
+ $ CALL DCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
*
-* Interchange rows KK and KP in last KK columns of A and W
+* Interchange rows KK and KP in last K+1 to N columns of A
+* (columns K (or K and K-1 for 2-by-2 pivot) of A will be
+* later overwritten). Interchange rows KK and KP
+* in last KKW to NB columns of W.
*
- CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
+ IF( K.LT.N )
+ $ CALL DSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ),
+ $ LDA )
CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
$ LDW )
END IF
*
IF( KSTEP.EQ.1 ) THEN
*
-* 1-by-1 pivot block D(k): column KW of W now holds
+* 1-by-1 pivot block D(k): column kw of W now holds
*
-* W(k) = U(k)*D(k)
+* W(kw) = U(k)*D(k),
*
* where U(k) is the k-th column of U
*
-* Store U(k) in column k of A
+* Store subdiag. elements of column U(k)
+* and 1-by-1 block D(k) in column k of A.
+* NOTE: Diagonal element U(k,k) is a UNIT element
+* and not stored.
+* A(k,k) := D(k,k) = W(k,kw)
+* A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
*
CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
R1 = ONE / A( K, K )
CALL DSCAL( K-1, R1, A( 1, K ), 1 )
+*
ELSE
*
-* 2-by-2 pivot block D(k): columns KW and KW-1 of W now
-* hold
+* 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
*
-* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
+* ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
*
* where U(k) and U(k-1) are the k-th and (k-1)-th columns
* of U
*
+* Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
+* block D(k-1:k,k-1:k) in columns k-1 and k of A.
+* NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
+* block and not stored.
+* A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
+* A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
+* = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
+*
IF( K.GT.2 ) THEN
*
-* Store U(k) and U(k-1) in columns k and k-1 of A
+* Compose the columns of the inverse of 2-by-2 pivot
+* block D in the following way to reduce the number
+* of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by
+* this inverse
+*
+* D**(-1) = ( d11 d21 )**(-1) =
+* ( d21 d22 )
+*
+* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
+* ( (-d21 ) ( d11 ) )
+*
+* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
+*
+* * ( ( d22/d21 ) ( -1 ) ) =
+* ( ( -1 ) ( d11/d21 ) )
+*
+* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
+* ( ( -1 ) ( D22 ) )
+*
+* = 1/d21 * T * ( ( D11 ) ( -1 ) )
+* ( ( -1 ) ( D22 ) )
+*
+* = D21 * ( ( D11 ) ( -1 ) )
+* ( ( -1 ) ( D22 ) )
*
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 ) )
@@ -370,7 +449,9 @@
A( K-1, K-1 ) = W( K-1, KW-1 )
A( K-1, K ) = W( K-1, KW )
A( K, K ) = W( K, KW )
+*
END IF
+*
END IF
*
* Store details of the interchanges in IPIV
@@ -414,20 +495,28 @@
50 CONTINUE
*
* Put U12 in standard form by partially undoing the interchanges
-* in columns k+1:n
+* in columns k+1:n looping backwards from k+1 to n
*
J = K + 1
60 CONTINUE
- JJ = J
- JP = IPIV( J )
- IF( JP.LT.0 ) THEN
- JP = -JP
+*
+* Undo the interchanges (if any) of rows JJ and JP at each
+* step J
+*
+* (Here, J is a diagonal index)
+ JJ = J
+ JP = IPIV( J )
+ IF( JP.LT.0 ) THEN
+ JP = -JP
+* (Here, J is a diagonal index)
+ J = J + 1
+ END IF
+* (NOTE: Here, J is used to determine row length. Length N-J+1
+* of the rows to swap back doesn't include diagonal element)
J = J + 1
- END IF
- J = J + 1
- IF( JP.NE.JJ .AND. J.LE.N )
- $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
- IF( J.LE.N )
+ IF( JP.NE.JJ .AND. J.LE.N )
+ $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
+ IF( J.LT.N )
$ GO TO 60
*
* Set KB to the number of columns factorized
@@ -464,7 +553,8 @@
ABSAKK = ABS( W( K, K ) )
*
* IMAX is the row-index of the largest off-diagonal element in
-* column K, and COLMAX is its absolute value
+* column K, and COLMAX is its absolute value.
+* Determine both COLMAX and IMAX.
*
IF( K.LT.N ) THEN
IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
@@ -475,7 +565,7 @@
*
IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
*
-* Column K is zero: set INFO and continue
+* Column K is zero or underflow: set INFO and continue
*
IF( INFO.EQ.0 )
$ INFO = K
@@ -518,7 +608,7 @@
*
KP = IMAX
*
-* copy column K+1 of W to column K
+* copy column K+1 of W to column K of W
*
CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
ELSE
@@ -531,21 +621,35 @@
END IF
END IF
*
+* ============================================================
+*
+* KK is the column of A where pivoting step stopped
+*
KK = K + KSTEP - 1
*
-* Updated column KP is already stored in column KK of W
+* Interchange rows and columns KP and KK.
+* Updated column KP is already stored in column KK of W.
*
IF( KP.NE.KK ) THEN
*
-* Copy non-updated column KK to column KP
+* Copy non-updated column KK to column KP of submatrix A
+* at step K. No need to copy element into column K
+* (or K and K+1 for 2-by-2 pivot) of A, since these columns
+* will be later overwritten.
*
- A( KP, K ) = A( KK, K )
- CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
- CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
+ A( KP, KP ) = A( KK, KK )
+ CALL DCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
+ $ LDA )
+ IF( KP.LT.N )
+ $ CALL DCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
*
-* Interchange rows KK and KP in first KK columns of A and W
+* Interchange rows KK and KP in first K-1 columns of A
+* (columns K (or K and K+1 for 2-by-2 pivot) of A will be
+* later overwritten). Interchange rows KK and KP
+* in first KK columns of W.
*
- CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
+ IF( K.GT.1 )
+ $ CALL DSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
END IF
*
@@ -553,17 +657,23 @@
*
* 1-by-1 pivot block D(k): column k of W now holds
*
-* W(k) = L(k)*D(k)
+* W(k) = L(k)*D(k),
*
* where L(k) is the k-th column of L
*
-* Store L(k) in column k of A
+* Store subdiag. elements of column L(k)
+* and 1-by-1 block D(k) in column k of A.
+* (NOTE: Diagonal element L(k,k) is a UNIT element
+* and not stored)
+* A(k,k) := D(k,k) = W(k,k)
+* A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
*
CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
IF( K.LT.N ) THEN
R1 = ONE / A( K, K )
CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
END IF
+*
ELSE
*
* 2-by-2 pivot block D(k): columns k and k+1 of W now hold
@@ -573,15 +683,51 @@
* where L(k) and L(k+1) are the k-th and (k+1)-th columns
* of L
*
+* Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
+* block D(k:k+1,k:k+1) in columns k and k+1 of A.
+* (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
+* block and not stored)
+* A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
+* A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
+* = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
+*
IF( K.LT.N-1 ) THEN
*
-* Store L(k) and L(k+1) in columns k and k+1 of A
+* Compose the columns of the inverse of 2-by-2 pivot
+* block D in the following way to reduce the number
+* of FLOPS when we myltiply panel ( W(k) W(k+1) ) by
+* this inverse
+*
+* D**(-1) = ( d11 d21 )**(-1) =
+* ( d21 d22 )
+*
+* = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
+* ( (-d21 ) ( d11 ) )
+*
+* = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
+*
+* * ( ( d22/d21 ) ( -1 ) ) =
+* ( ( -1 ) ( d11/d21 ) )
+*
+* = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) =
+* ( ( -1 ) ( D22 ) )
+*
+* = 1/d21 * T * ( ( D11 ) ( -1 ) )
+* ( ( -1 ) ( D22 ) )
+*
+* = D21 * ( ( D11 ) ( -1 ) )
+* ( ( -1 ) ( D22 ) )
*
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 ) )
@@ -593,7 +739,9 @@
A( K, K ) = W( K, K )
A( K+1, K ) = W( K+1, K )
A( K+1, K+1 ) = W( K+1, K+1 )
+*
END IF
+*
END IF
*
* Store details of the interchanges in IPIV
@@ -638,20 +786,28 @@
110 CONTINUE
*
* Put L21 in standard form by partially undoing the interchanges
-* in columns 1:k-1
+* of rows in columns 1:k-1 looping backwards from k-1 to 1
*
J = K - 1
120 CONTINUE
- JJ = J
- JP = IPIV( J )
- IF( JP.LT.0 ) THEN
- JP = -JP
+*
+* Undo the interchanges (if any) of rows JJ and JP at each
+* step J
+*
+* (Here, J is a diagonal index)
+ JJ = J
+ JP = IPIV( J )
+ IF( JP.LT.0 ) THEN
+ JP = -JP
+* (Here, J is a diagonal index)
+ J = J - 1
+ END IF
+* (NOTE: Here, J is used to determine row length. Length J
+* of the rows to swap back doesn't include diagonal element)
J = J - 1
- END IF
- J = J - 1
- IF( JP.NE.JJ .AND. J.GE.1 )
- $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
- IF( J.GE.1 )
+ IF( JP.NE.JJ .AND. J.GE.1 )
+ $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
+ IF( J.GT.1 )
$ GO TO 120
*
* Set KB to the number of columns factorized