--- rpl/lapack/lapack/zlaqps.f 2010/04/21 13:45:34 1.2 +++ rpl/lapack/lapack/zlaqps.f 2011/11/21 20:43:16 1.10 @@ -1,10 +1,186 @@ +*> \brief \b ZLAQPS +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLAQPS + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, +* VN2, AUXV, F, LDF ) +* +* .. Scalar Arguments .. +* INTEGER KB, LDA, LDF, M, N, NB, OFFSET +* .. +* .. Array Arguments .. +* INTEGER JPVT( * ) +* DOUBLE PRECISION VN1( * ), VN2( * ) +* COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLAQPS computes a step of QR factorization with column pivoting +*> of a complex M-by-N matrix A by using Blas-3. It tries to factorize +*> NB columns from A starting from the row OFFSET+1, and updates all +*> of the matrix with Blas-3 xGEMM. +*> +*> In some cases, due to catastrophic cancellations, it cannot +*> factorize NB columns. Hence, the actual number of factorized +*> columns is returned in KB. +*> +*> Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix A. N >= 0 +*> \endverbatim +*> +*> \param[in] OFFSET +*> \verbatim +*> OFFSET is INTEGER +*> The number of rows of A that have been factorized in +*> previous steps. +*> \endverbatim +*> +*> \param[in] NB +*> \verbatim +*> NB is INTEGER +*> The number of columns to factorize. +*> \endverbatim +*> +*> \param[out] KB +*> \verbatim +*> KB is INTEGER +*> The number of columns actually factorized. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, block A(OFFSET+1:M,1:KB) is the triangular +*> factor obtained and block A(1:OFFSET,1:N) has been +*> accordingly pivoted, but no factorized. +*> The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has +*> been updated. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] JPVT +*> \verbatim +*> JPVT is INTEGER array, dimension (N) +*> JPVT(I) = K <==> Column K of the full matrix A has been +*> permuted into position I in AP. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX*16 array, dimension (KB) +*> The scalar factors of the elementary reflectors. +*> \endverbatim +*> +*> \param[in,out] VN1 +*> \verbatim +*> VN1 is DOUBLE PRECISION array, dimension (N) +*> The vector with the partial column norms. +*> \endverbatim +*> +*> \param[in,out] VN2 +*> \verbatim +*> VN2 is DOUBLE PRECISION array, dimension (N) +*> The vector with the exact column norms. +*> \endverbatim +*> +*> \param[in,out] AUXV +*> \verbatim +*> AUXV is COMPLEX*16 array, dimension (NB) +*> Auxiliar vector. +*> \endverbatim +*> +*> \param[in,out] F +*> \verbatim +*> F is COMPLEX*16 array, dimension (LDF,NB) +*> Matrix F**H = L * Y**H * A. +*> \endverbatim +*> +*> \param[in] LDF +*> \verbatim +*> LDF is INTEGER +*> The leading dimension of the array F. LDF >= max(1,N). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERauxiliary +* +*> \par Contributors: +* ================== +*> +*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain +*> X. Sun, Computer Science Dept., Duke University, USA +*> \n +*> Partial column norm updating strategy modified on April 2011 +*> Z. Drmac and Z. Bujanovic, Dept. of Mathematics, +*> University of Zagreb, Croatia. +* +*> \par References: +* ================ +*> +*> LAPACK Working Note 176 +* +*> \htmlonly +*> [PDF] +*> \endhtmlonly +* +* ===================================================================== SUBROUTINE ZLAQPS( M, N, OFFSET, NB, KB, A, LDA, JPVT, TAU, VN1, $ VN2, AUXV, F, LDF ) * -* -- LAPACK auxiliary routine (version 3.2) -- +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER KB, LDA, LDF, M, N, NB, OFFSET @@ -15,79 +191,6 @@ COMPLEX*16 A( LDA, * ), AUXV( * ), F( LDF, * ), TAU( * ) * .. * -* Purpose -* ======= -* -* ZLAQPS computes a step of QR factorization with column pivoting -* of a complex M-by-N matrix A by using Blas-3. It tries to factorize -* NB columns from A starting from the row OFFSET+1, and updates all -* of the matrix with Blas-3 xGEMM. -* -* In some cases, due to catastrophic cancellations, it cannot -* factorize NB columns. Hence, the actual number of factorized -* columns is returned in KB. -* -* Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. -* -* Arguments -* ========= -* -* M (input) INTEGER -* The number of rows of the matrix A. M >= 0. -* -* N (input) INTEGER -* The number of columns of the matrix A. N >= 0 -* -* OFFSET (input) INTEGER -* The number of rows of A that have been factorized in -* previous steps. -* -* NB (input) INTEGER -* The number of columns to factorize. -* -* KB (output) INTEGER -* The number of columns actually factorized. -* -* A (input/output) COMPLEX*16 array, dimension (LDA,N) -* On entry, the M-by-N matrix A. -* On exit, block A(OFFSET+1:M,1:KB) is the triangular -* factor obtained and block A(1:OFFSET,1:N) has been -* accordingly pivoted, but no factorized. -* The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has -* been updated. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,M). -* -* JPVT (input/output) INTEGER array, dimension (N) -* JPVT(I) = K <==> Column K of the full matrix A has been -* permuted into position I in AP. -* -* TAU (output) COMPLEX*16 array, dimension (KB) -* The scalar factors of the elementary reflectors. -* -* VN1 (input/output) DOUBLE PRECISION array, dimension (N) -* The vector with the partial column norms. -* -* VN2 (input/output) DOUBLE PRECISION array, dimension (N) -* The vector with the exact column norms. -* -* AUXV (input/output) COMPLEX*16 array, dimension (NB) -* Auxiliar vector. -* -* F (input/output) COMPLEX*16 array, dimension (LDF,NB) -* Matrix F' = L*Y'*A. -* -* LDF (input) INTEGER -* The leading dimension of the array F. LDF >= max(1,N). -* -* Further Details -* =============== -* -* Based on contributions by -* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain -* X. Sun, Computer Science Dept., Duke University, USA -* * ===================================================================== * * .. Parameters .. @@ -103,7 +206,7 @@ COMPLEX*16 AKK * .. * .. External Subroutines .. - EXTERNAL ZGEMM, ZGEMV, ZLARFP, ZSWAP + EXTERNAL ZGEMM, ZGEMV, ZLARFG, ZSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, NINT, SQRT @@ -141,7 +244,7 @@ END IF * * Apply previous Householder reflectors to column K: -* A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. +* A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)**H. * IF( K.GT.1 ) THEN DO 20 J = 1, K - 1 @@ -157,9 +260,9 @@ * Generate elementary reflector H(k). * IF( RK.LT.M ) THEN - CALL ZLARFP( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) + CALL ZLARFG( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) ) ELSE - CALL ZLARFP( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) + CALL ZLARFG( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) ) END IF * AKK = A( RK, K ) @@ -167,7 +270,7 @@ * * Compute Kth column of F: * -* Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). +* Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K). * IF( K.LT.N ) THEN CALL ZGEMV( 'Conjugate transpose', M-RK+1, N-K, TAU( K ), @@ -182,7 +285,7 @@ 40 CONTINUE * * Incremental updating of F: -* F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' +* F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)**H * *A(RK:M,K). * IF( K.GT.1 ) THEN @@ -195,7 +298,7 @@ END IF * * Update the current row of A: -* A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. +* A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)**H. * IF( K.LT.N ) THEN CALL ZGEMM( 'No transpose', 'Conjugate transpose', 1, N-K, @@ -236,7 +339,7 @@ * * Apply the block reflector to the rest of the matrix: * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - -* A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. +* A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)**H. * IF( KB.LT.MIN( N, M-OFFSET ) ) THEN CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-RK, N-KB,