1*> \brief \b ZQLT02 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZQLT02( M, N, K, A, AF, Q, L, LDA, TAU, WORK, LWORK, 12* RWORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* INTEGER K, LDA, LWORK, M, N 16* .. 17* .. Array Arguments .. 18* DOUBLE PRECISION RESULT( * ), RWORK( * ) 19* COMPLEX*16 A( LDA, * ), AF( LDA, * ), L( LDA, * ), 20* $ Q( LDA, * ), TAU( * ), WORK( LWORK ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZQLT02 tests ZUNGQL, which generates an m-by-n matrix Q with 30*> orthonornmal columns that is defined as the product of k elementary 31*> reflectors. 32*> 33*> Given the QL factorization of an m-by-n matrix A, ZQLT02 generates 34*> the orthogonal matrix Q defined by the factorization of the last k 35*> columns of A; it compares L(m-n+1:m,n-k+1:n) with 36*> Q(1:m,m-n+1:m)'*A(1:m,n-k+1:n), and checks that the columns of Q are 37*> orthonormal. 38*> \endverbatim 39* 40* Arguments: 41* ========== 42* 43*> \param[in] M 44*> \verbatim 45*> M is INTEGER 46*> The number of rows of the matrix Q to be generated. M >= 0. 47*> \endverbatim 48*> 49*> \param[in] N 50*> \verbatim 51*> N is INTEGER 52*> The number of columns of the matrix Q to be generated. 53*> M >= N >= 0. 54*> \endverbatim 55*> 56*> \param[in] K 57*> \verbatim 58*> K is INTEGER 59*> The number of elementary reflectors whose product defines the 60*> matrix Q. N >= K >= 0. 61*> \endverbatim 62*> 63*> \param[in] A 64*> \verbatim 65*> A is COMPLEX*16 array, dimension (LDA,N) 66*> The m-by-n matrix A which was factorized by ZQLT01. 67*> \endverbatim 68*> 69*> \param[in] AF 70*> \verbatim 71*> AF is COMPLEX*16 array, dimension (LDA,N) 72*> Details of the QL factorization of A, as returned by ZGEQLF. 73*> See ZGEQLF for further details. 74*> \endverbatim 75*> 76*> \param[out] Q 77*> \verbatim 78*> Q is COMPLEX*16 array, dimension (LDA,N) 79*> \endverbatim 80*> 81*> \param[out] L 82*> \verbatim 83*> L is COMPLEX*16 array, dimension (LDA,N) 84*> \endverbatim 85*> 86*> \param[in] LDA 87*> \verbatim 88*> LDA is INTEGER 89*> The leading dimension of the arrays A, AF, Q and L. LDA >= M. 90*> \endverbatim 91*> 92*> \param[in] TAU 93*> \verbatim 94*> TAU is COMPLEX*16 array, dimension (N) 95*> The scalar factors of the elementary reflectors corresponding 96*> to the QL factorization in AF. 97*> \endverbatim 98*> 99*> \param[out] WORK 100*> \verbatim 101*> WORK is COMPLEX*16 array, dimension (LWORK) 102*> \endverbatim 103*> 104*> \param[in] LWORK 105*> \verbatim 106*> LWORK is INTEGER 107*> The dimension of the array WORK. 108*> \endverbatim 109*> 110*> \param[out] RWORK 111*> \verbatim 112*> RWORK is DOUBLE PRECISION array, dimension (M) 113*> \endverbatim 114*> 115*> \param[out] RESULT 116*> \verbatim 117*> RESULT is DOUBLE PRECISION array, dimension (2) 118*> The test ratios: 119*> RESULT(1) = norm( L - Q'*A ) / ( M * norm(A) * EPS ) 120*> RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) 121*> \endverbatim 122* 123* Authors: 124* ======== 125* 126*> \author Univ. of Tennessee 127*> \author Univ. of California Berkeley 128*> \author Univ. of Colorado Denver 129*> \author NAG Ltd. 130* 131*> \date November 2011 132* 133*> \ingroup complex16_lin 134* 135* ===================================================================== 136 SUBROUTINE ZQLT02( M, N, K, A, AF, Q, L, LDA, TAU, WORK, LWORK, 137 $ RWORK, RESULT ) 138* 139* -- LAPACK test routine (version 3.4.0) -- 140* -- LAPACK is a software package provided by Univ. of Tennessee, -- 141* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 142* November 2011 143* 144* .. Scalar Arguments .. 145 INTEGER K, LDA, LWORK, M, N 146* .. 147* .. Array Arguments .. 148 DOUBLE PRECISION RESULT( * ), RWORK( * ) 149 COMPLEX*16 A( LDA, * ), AF( LDA, * ), L( LDA, * ), 150 $ Q( LDA, * ), TAU( * ), WORK( LWORK ) 151* .. 152* 153* ===================================================================== 154* 155* .. Parameters .. 156 DOUBLE PRECISION ZERO, ONE 157 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 158 COMPLEX*16 ROGUE 159 PARAMETER ( ROGUE = ( -1.0D+10, -1.0D+10 ) ) 160* .. 161* .. Local Scalars .. 162 INTEGER INFO 163 DOUBLE PRECISION ANORM, EPS, RESID 164* .. 165* .. External Functions .. 166 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY 167 EXTERNAL DLAMCH, ZLANGE, ZLANSY 168* .. 169* .. External Subroutines .. 170 EXTERNAL ZGEMM, ZHERK, ZLACPY, ZLASET, ZUNGQL 171* .. 172* .. Intrinsic Functions .. 173 INTRINSIC DBLE, DCMPLX, MAX 174* .. 175* .. Scalars in Common .. 176 CHARACTER*32 SRNAMT 177* .. 178* .. Common blocks .. 179 COMMON / SRNAMC / SRNAMT 180* .. 181* .. Executable Statements .. 182* 183* Quick return if possible 184* 185 IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN 186 RESULT( 1 ) = ZERO 187 RESULT( 2 ) = ZERO 188 RETURN 189 END IF 190* 191 EPS = DLAMCH( 'Epsilon' ) 192* 193* Copy the last k columns of the factorization to the array Q 194* 195 CALL ZLASET( 'Full', M, N, ROGUE, ROGUE, Q, LDA ) 196 IF( K.LT.M ) 197 $ CALL ZLACPY( 'Full', M-K, K, AF( 1, N-K+1 ), LDA, 198 $ Q( 1, N-K+1 ), LDA ) 199 IF( K.GT.1 ) 200 $ CALL ZLACPY( 'Upper', K-1, K-1, AF( M-K+1, N-K+2 ), LDA, 201 $ Q( M-K+1, N-K+2 ), LDA ) 202* 203* Generate the last n columns of the matrix Q 204* 205 SRNAMT = 'ZUNGQL' 206 CALL ZUNGQL( M, N, K, Q, LDA, TAU( N-K+1 ), WORK, LWORK, INFO ) 207* 208* Copy L(m-n+1:m,n-k+1:n) 209* 210 CALL ZLASET( 'Full', N, K, DCMPLX( ZERO ), DCMPLX( ZERO ), 211 $ L( M-N+1, N-K+1 ), LDA ) 212 CALL ZLACPY( 'Lower', K, K, AF( M-K+1, N-K+1 ), LDA, 213 $ L( M-K+1, N-K+1 ), LDA ) 214* 215* Compute L(m-n+1:m,n-k+1:n) - Q(1:m,m-n+1:m)' * A(1:m,n-k+1:n) 216* 217 CALL ZGEMM( 'Conjugate transpose', 'No transpose', N, K, M, 218 $ DCMPLX( -ONE ), Q, LDA, A( 1, N-K+1 ), LDA, 219 $ DCMPLX( ONE ), L( M-N+1, N-K+1 ), LDA ) 220* 221* Compute norm( L - Q'*A ) / ( M * norm(A) * EPS ) . 222* 223 ANORM = ZLANGE( '1', M, K, A( 1, N-K+1 ), LDA, RWORK ) 224 RESID = ZLANGE( '1', N, K, L( M-N+1, N-K+1 ), LDA, RWORK ) 225 IF( ANORM.GT.ZERO ) THEN 226 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS 227 ELSE 228 RESULT( 1 ) = ZERO 229 END IF 230* 231* Compute I - Q'*Q 232* 233 CALL ZLASET( 'Full', N, N, DCMPLX( ZERO ), DCMPLX( ONE ), L, LDA ) 234 CALL ZHERK( 'Upper', 'Conjugate transpose', N, M, -ONE, Q, LDA, 235 $ ONE, L, LDA ) 236* 237* Compute norm( I - Q'*Q ) / ( M * EPS ) . 238* 239 RESID = ZLANSY( '1', 'Upper', N, L, LDA, RWORK ) 240* 241 RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS 242* 243 RETURN 244* 245* End of ZQLT02 246* 247 END 248