1*> \brief \b DQLT02 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 DQLT02( 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 A( LDA, * ), AF( LDA, * ), L( LDA, * ), 19* $ Q( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), 20* $ WORK( LWORK ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> DQLT02 tests DORGQL, 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, DQLT02 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 DOUBLE PRECISION array, dimension (LDA,N) 66*> The m-by-n matrix A which was factorized by DQLT01. 67*> \endverbatim 68*> 69*> \param[in] AF 70*> \verbatim 71*> AF is DOUBLE PRECISION array, dimension (LDA,N) 72*> Details of the QL factorization of A, as returned by DGEQLF. 73*> See DGEQLF for further details. 74*> \endverbatim 75*> 76*> \param[out] Q 77*> \verbatim 78*> Q is DOUBLE PRECISION array, dimension (LDA,N) 79*> \endverbatim 80*> 81*> \param[out] L 82*> \verbatim 83*> L is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 December 2016 132* 133*> \ingroup double_lin 134* 135* ===================================================================== 136 SUBROUTINE DQLT02( M, N, K, A, AF, Q, L, LDA, TAU, WORK, LWORK, 137 $ RWORK, RESULT ) 138* 139* -- LAPACK test routine (version 3.7.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* December 2016 143* 144* .. Scalar Arguments .. 145 INTEGER K, LDA, LWORK, M, N 146* .. 147* .. Array Arguments .. 148 DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), L( LDA, * ), 149 $ Q( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), 150 $ WORK( LWORK ) 151* .. 152* 153* ===================================================================== 154* 155* .. Parameters .. 156 DOUBLE PRECISION ZERO, ONE 157 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 158 DOUBLE PRECISION ROGUE 159 PARAMETER ( ROGUE = -1.0D+10 ) 160* .. 161* .. Local Scalars .. 162 INTEGER INFO 163 DOUBLE PRECISION ANORM, EPS, RESID 164* .. 165* .. External Functions .. 166 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 167 EXTERNAL DLAMCH, DLANGE, DLANSY 168* .. 169* .. External Subroutines .. 170 EXTERNAL DGEMM, DLACPY, DLASET, DORGQL, DSYRK 171* .. 172* .. Intrinsic Functions .. 173 INTRINSIC DBLE, 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 DLASET( 'Full', M, N, ROGUE, ROGUE, Q, LDA ) 196 IF( K.LT.M ) 197 $ CALL DLACPY( 'Full', M-K, K, AF( 1, N-K+1 ), LDA, 198 $ Q( 1, N-K+1 ), LDA ) 199 IF( K.GT.1 ) 200 $ CALL DLACPY( '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 = 'DORGQL' 206 CALL DORGQL( 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 DLASET( 'Full', N, K, ZERO, ZERO, L( M-N+1, N-K+1 ), LDA ) 211 CALL DLACPY( 'Lower', K, K, AF( M-K+1, N-K+1 ), LDA, 212 $ L( M-K+1, N-K+1 ), LDA ) 213* 214* Compute L(m-n+1:m,n-k+1:n) - Q(1:m,m-n+1:m)' * A(1:m,n-k+1:n) 215* 216 CALL DGEMM( 'Transpose', 'No transpose', N, K, M, -ONE, Q, LDA, 217 $ A( 1, N-K+1 ), LDA, ONE, L( M-N+1, N-K+1 ), LDA ) 218* 219* Compute norm( L - Q'*A ) / ( M * norm(A) * EPS ) . 220* 221 ANORM = DLANGE( '1', M, K, A( 1, N-K+1 ), LDA, RWORK ) 222 RESID = DLANGE( '1', N, K, L( M-N+1, N-K+1 ), LDA, RWORK ) 223 IF( ANORM.GT.ZERO ) THEN 224 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M ) ) ) / ANORM ) / EPS 225 ELSE 226 RESULT( 1 ) = ZERO 227 END IF 228* 229* Compute I - Q'*Q 230* 231 CALL DLASET( 'Full', N, N, ZERO, ONE, L, LDA ) 232 CALL DSYRK( 'Upper', 'Transpose', N, M, -ONE, Q, LDA, ONE, L, 233 $ LDA ) 234* 235* Compute norm( I - Q'*Q ) / ( M * EPS ) . 236* 237 RESID = DLANSY( '1', 'Upper', N, L, LDA, RWORK ) 238* 239 RESULT( 2 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / EPS 240* 241 RETURN 242* 243* End of DQLT02 244* 245 END 246