--- rpl/lapack/lapack/zlahr2.f 2010/08/13 21:04:08 1.6
+++ rpl/lapack/lapack/zlahr2.f 2023/08/07 08:39:29 1.19
@@ -1,9 +1,187 @@
+*> \brief \b ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZLAHR2 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
+*
+* .. Scalar Arguments ..
+* INTEGER K, LDA, LDT, LDY, N, NB
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), T( LDT, NB ), TAU( NB ),
+* $ Y( LDY, NB )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1)
+*> matrix A so that elements below the k-th subdiagonal are zero. The
+*> reduction is performed by an unitary similarity transformation
+*> Q**H * A * Q. The routine returns the matrices V and T which determine
+*> Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T.
+*>
+*> This is an auxiliary routine called by ZGEHRD.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*> K is INTEGER
+*> The offset for the reduction. Elements below the k-th
+*> subdiagonal in the first NB columns are reduced to zero.
+*> K < N.
+*> \endverbatim
+*>
+*> \param[in] NB
+*> \verbatim
+*> NB is INTEGER
+*> The number of columns to be reduced.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N-K+1)
+*> On entry, the n-by-(n-k+1) general matrix A.
+*> On exit, the elements on and above the k-th subdiagonal in
+*> the first NB columns are overwritten with the corresponding
+*> elements of the reduced matrix; the elements below the k-th
+*> subdiagonal, with the array TAU, represent the matrix Q as a
+*> product of elementary reflectors. The other columns of A are
+*> unchanged. See Further Details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*> TAU is COMPLEX*16 array, dimension (NB)
+*> The scalar factors of the elementary reflectors. See Further
+*> Details.
+*> \endverbatim
+*>
+*> \param[out] T
+*> \verbatim
+*> T is COMPLEX*16 array, dimension (LDT,NB)
+*> The upper triangular matrix T.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*> LDT is INTEGER
+*> The leading dimension of the array T. LDT >= NB.
+*> \endverbatim
+*>
+*> \param[out] Y
+*> \verbatim
+*> Y is COMPLEX*16 array, dimension (LDY,NB)
+*> The n-by-nb matrix Y.
+*> \endverbatim
+*>
+*> \param[in] LDY
+*> \verbatim
+*> LDY is INTEGER
+*> The leading dimension of the array Y. LDY >= N.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> The matrix Q is represented as a product of nb elementary reflectors
+*>
+*> Q = H(1) H(2) . . . H(nb).
+*>
+*> Each H(i) has the form
+*>
+*> H(i) = I - tau * v * v**H
+*>
+*> where tau is a complex scalar, and v is a complex vector with
+*> v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
+*> A(i+k+1:n,i), and tau in TAU(i).
+*>
+*> The elements of the vectors v together form the (n-k+1)-by-nb matrix
+*> V which is needed, with T and Y, to apply the transformation to the
+*> unreduced part of the matrix, using an update of the form:
+*> A := (I - V*T*V**H) * (A - Y*V**H).
+*>
+*> The contents of A on exit are illustrated by the following example
+*> with n = 7, k = 3 and nb = 2:
+*>
+*> ( a a a a a )
+*> ( a a a a a )
+*> ( a a a a a )
+*> ( h h a a a )
+*> ( v1 h a a a )
+*> ( v1 v2 a a a )
+*> ( v1 v2 a a a )
+*>
+*> where a denotes an element of the original matrix A, h denotes a
+*> modified element of the upper Hessenberg matrix H, and vi denotes an
+*> element of the vector defining H(i).
+*>
+*> This subroutine is a slight modification of LAPACK-3.0's ZLAHRD
+*> incorporating improvements proposed by Quintana-Orti and Van de
+*> Gejin. Note that the entries of A(1:K,2:NB) differ from those
+*> returned by the original LAPACK-3.0's ZLAHRD routine. (This
+*> subroutine is not backward compatible with LAPACK-3.0's ZLAHRD.)
+*> \endverbatim
+*
+*> \par References:
+* ================
+*>
+*> Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
+*> performance of reduction to Hessenberg form," ACM Transactions on
+*> Mathematical Software, 32(2):180-194, June 2006.
+*>
+* =====================================================================
SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
*
-* -- LAPACK auxiliary routine (version 3.2.1) --
+* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* -- April 2009 --
*
* .. Scalar Arguments ..
INTEGER K, LDA, LDT, LDY, N, NB
@@ -13,112 +191,11 @@
$ Y( LDY, NB )
* ..
*
-* Purpose
-* =======
-*
-* ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1)
-* matrix A so that elements below the k-th subdiagonal are zero. The
-* reduction is performed by an unitary similarity transformation
-* Q' * A * Q. The routine returns the matrices V and T which determine
-* Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
-*
-* This is an auxiliary routine called by ZGEHRD.
-*
-* Arguments
-* =========
-*
-* N (input) INTEGER
-* The order of the matrix A.
-*
-* K (input) INTEGER
-* The offset for the reduction. Elements below the k-th
-* subdiagonal in the first NB columns are reduced to zero.
-* K < N.
-*
-* NB (input) INTEGER
-* The number of columns to be reduced.
-*
-* A (input/output) COMPLEX*16 array, dimension (LDA,N-K+1)
-* On entry, the n-by-(n-k+1) general matrix A.
-* On exit, the elements on and above the k-th subdiagonal in
-* the first NB columns are overwritten with the corresponding
-* elements of the reduced matrix; the elements below the k-th
-* subdiagonal, with the array TAU, represent the matrix Q as a
-* product of elementary reflectors. The other columns of A are
-* unchanged. See Further Details.
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,N).
-*
-* TAU (output) COMPLEX*16 array, dimension (NB)
-* The scalar factors of the elementary reflectors. See Further
-* Details.
-*
-* T (output) COMPLEX*16 array, dimension (LDT,NB)
-* The upper triangular matrix T.
-*
-* LDT (input) INTEGER
-* The leading dimension of the array T. LDT >= NB.
-*
-* Y (output) COMPLEX*16 array, dimension (LDY,NB)
-* The n-by-nb matrix Y.
-*
-* LDY (input) INTEGER
-* The leading dimension of the array Y. LDY >= N.
-*
-* Further Details
-* ===============
-*
-* The matrix Q is represented as a product of nb elementary reflectors
-*
-* Q = H(1) H(2) . . . H(nb).
-*
-* Each H(i) has the form
-*
-* H(i) = I - tau * v * v'
-*
-* where tau is a complex scalar, and v is a complex vector with
-* v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
-* A(i+k+1:n,i), and tau in TAU(i).
-*
-* The elements of the vectors v together form the (n-k+1)-by-nb matrix
-* V which is needed, with T and Y, to apply the transformation to the
-* unreduced part of the matrix, using an update of the form:
-* A := (I - V*T*V') * (A - Y*V').
-*
-* The contents of A on exit are illustrated by the following example
-* with n = 7, k = 3 and nb = 2:
-*
-* ( a a a a a )
-* ( a a a a a )
-* ( a a a a a )
-* ( h h a a a )
-* ( v1 h a a a )
-* ( v1 v2 a a a )
-* ( v1 v2 a a a )
-*
-* where a denotes an element of the original matrix A, h denotes a
-* modified element of the upper Hessenberg matrix H, and vi denotes an
-* element of the vector defining H(i).
-*
-* This subroutine is a slight modification of LAPACK-3.0's DLAHRD
-* incorporating improvements proposed by Quintana-Orti and Van de
-* Gejin. Note that the entries of A(1:K,2:NB) differ from those
-* returned by the original LAPACK-3.0's DLAHRD routine. (This
-* subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
-*
-* References
-* ==========
-*
-* Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
-* performance of reduction to Hessenberg form," ACM Transactions on
-* Mathematical Software, 32(2):180-194, June 2006.
-*
* =====================================================================
*
* .. Parameters ..
COMPLEX*16 ZERO, ONE
- PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
+ PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
$ ONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
@@ -144,14 +221,14 @@
*
* Update A(K+1:N,I)
*
-* Update I-th column of A - Y * V'
+* Update I-th column of A - Y * V**H
*
- CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
+ CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
$ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
- CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
+ CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
*
-* Apply I - V * T' * V' to this column (call it b) from the
+* Apply I - V * T**H * V**H to this column (call it b) from the
* left, using the last column of T as workspace
*
* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
@@ -159,34 +236,34 @@
*
* where V1 is unit lower triangular
*
-* w := V1' * b1
+* w := V1**H * b1
*
CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
- CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT',
+ CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT',
$ I-1, A( K+1, 1 ),
$ LDA, T( 1, NB ), 1 )
*
-* w := w + V2'*b2
+* w := w + V2**H * b2
*
- CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
+ CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
$ ONE, A( K+I, 1 ),
$ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
*
-* w := T'*w
+* w := T**H * w
*
- CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT',
+ CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT',
$ I-1, T, LDT,
$ T( 1, NB ), 1 )
*
* b2 := b2 - V2*w
*
- CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE,
+ CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE,
$ A( K+I, 1 ),
$ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
*
* b1 := b1 - V1*w
*
- CALL ZTRMV( 'Lower', 'NO TRANSPOSE',
+ CALL ZTRMV( 'Lower', 'NO TRANSPOSE',
$ 'UNIT', I-1,
$ A( K+1, 1 ), LDA, T( 1, NB ), 1 )
CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
@@ -204,13 +281,13 @@
*
* Compute Y(K+1:N,I)
*
- CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1,
+ CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1,
$ ONE, A( K+1, I+1 ),
$ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
- CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
+ CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
$ ONE, A( K+I, 1 ), LDA,
$ A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
- CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE,
+ CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE,
$ Y( K+1, 1 ), LDY,
$ T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
CALL ZSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
@@ -218,7 +295,7 @@
* Compute T(1:I,I)
*
CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
- CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT',
+ CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT',
$ I-1, T, LDT,
$ T( 1, I ), 1 )
T( I, I ) = TAU( I )
@@ -229,15 +306,15 @@
* Compute Y(1:K,1:NB)
*
CALL ZLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
- CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE',
+ CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE',
$ 'UNIT', K, NB,
$ ONE, A( K+1, 1 ), LDA, Y, LDY )
IF( N.GT.K+NB )
- $ CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K,
+ $ CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K,
$ NB, N-K-NB, ONE,
$ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
$ LDY )
- CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE',
+ CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE',
$ 'NON-UNIT', K, NB,
$ ONE, T, LDT, Y, LDY )
*