From 62a53eb896fb85ffdf8760171b7868a2132a0188 Mon Sep 17 00:00:00 2001 From: igor175 Date: Thu, 8 Nov 2012 21:53:21 +0000 Subject: (s,d,c,z)lasyf.f and (c,z)lahef.f: added comments to the part where a column or 2 columns are updated at each step --- SRC/clahef.f | 132 ++++++++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 118 insertions(+), 14 deletions(-) (limited to 'SRC/clahef.f') diff --git a/SRC/clahef.f b/SRC/clahef.f index b0ffc35f..fc2c3392 100644 --- a/SRC/clahef.f +++ b/SRC/clahef.f @@ -333,8 +333,13 @@ KSTEP = 2 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 * * Interchange rows and columns KP and KK. @@ -368,40 +373,86 @@ * 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 +* (1) 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 CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) R1 = ONE / REAL( A( K, K ) ) CALL CSSCAL( K-1, R1, A( 1, K ), 1 ) * -* Conjugate W(k) +* (2) Conjugate column W(kw) * CALL CLACGV( K-1, W( 1, KW ), 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 +* +* (1) 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 cj(d21) )**(-1) = +* ( d21 d22 ) +* +* = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) = +* ( (-d21) ( d11 ) ) +* +* = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) * +* +* * ( d21*( d22/d21 ) conj(d21)( - 1 ) ) = +* ( ( -1 ) ( d11/conj(d21) ) ) +* +* = 1/(|d21|**2) * 1/(D22*D11-1) * +* +* * ( d21*( D11 ) conj(d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = ( conj(D21)*( D11 ) D21*( -1 ) ) +* ( ( -1 ) ( D22 ) ) * D21 = W( K-1, KW ) 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 )* @@ -415,11 +466,13 @@ A( K-1, K ) = W( K-1, KW ) A( K, K ) = W( K, KW ) * -* Conjugate W(k) and W(k-1) +* (2) Conjugate columns W(kw) and W(kw-1) * CALL CLACGV( K-1, W( 1, KW ), 1 ) CALL CLACGV( K-2, W( 1, KW-1 ), 1 ) +* END IF +* END IF * * Store details of the interchanges in IPIV @@ -592,6 +645,8 @@ KSTEP = 2 END IF END IF +* +* KK is the column of A where pivoting step stopped * KK = K + KSTEP - 1 * @@ -626,21 +681,27 @@ * * 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 +* (1) 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 CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) IF( K.LT.N ) THEN R1 = ONE / REAL( A( K, K ) ) CALL CSSCAL( N-K, R1, A( K+1, K ), 1 ) * -* Conjugate W(k) +* (2) Conjugate column W(k) * CALL CLACGV( N-K, W( K+1, K ), 1 ) END IF +* ELSE * * 2-by-2 pivot block D(k): columns k and k+1 of W now hold @@ -649,16 +710,57 @@ * * where L(k) and L(k+1) are the k-th and (k+1)-th columns * of L +* +* (1) 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(kw-1) W(kw) ) by +* this inverse +* +* D**(-1) = ( d11 cj(d21) )**(-1) = +* ( d21 d22 ) +* +* = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) = +* ( (-d21) ( d11 ) ) +* +* = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) * +* +* * ( d21*( d22/d21 ) conj(d21)( - 1 ) ) = +* ( ( -1 ) ( d11/conj(d21) ) ) +* +* = 1/(|d21|**2) * 1/(D22*D11-1) * +* +* * ( d21*( D11 ) conj(d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) = +* ( ( -1 ) ( D22 ) ) +* +* = ( conj(D21)*( D11 ) D21*( -1 ) ) +* ( ( -1 ) ( D22 ) ) * D21 = W( K+1, K ) 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 ) ) @@ -672,11 +774,13 @@ A( K+1, K ) = W( K+1, K ) A( K+1, K+1 ) = W( K+1, K+1 ) * -* Conjugate W(k) and W(k+1) +* (2) Conjugate columns W(k) and W(k+1) * CALL CLACGV( N-K, W( K+1, K ), 1 ) CALL CLACGV( N-K-1, W( K+2, K+1 ), 1 ) +* END IF +* END IF * * Store details of the interchanges in IPIV -- cgit v1.2.3