![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGRQTS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE ZGRQTS( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, 00012 * BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDB, LWORK, M, N, P 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION RESULT( 4 ), RWORK( * ) 00019 * COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00020 * $ BF( LDB, * ), BWK( LDB, * ), Q( LDA, * ), 00021 * $ R( LDA, * ), T( LDB, * ), TAUA( * ), TAUB( * ), 00022 * $ WORK( LWORK ), Z( LDB, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> ZGRQTS tests ZGGRQF, which computes the GRQ factorization of an 00032 *> M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] M 00039 *> \verbatim 00040 *> M is INTEGER 00041 *> The number of rows of the matrix A. M >= 0. 00042 *> \endverbatim 00043 *> 00044 *> \param[in] P 00045 *> \verbatim 00046 *> P is INTEGER 00047 *> The number of rows of the matrix B. P >= 0. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] N 00051 *> \verbatim 00052 *> N is INTEGER 00053 *> The number of columns of the matrices A and B. N >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] A 00057 *> \verbatim 00058 *> A is COMPLEX*16 array, dimension (LDA,N) 00059 *> The M-by-N matrix A. 00060 *> \endverbatim 00061 *> 00062 *> \param[out] AF 00063 *> \verbatim 00064 *> AF is COMPLEX*16 array, dimension (LDA,N) 00065 *> Details of the GRQ factorization of A and B, as returned 00066 *> by ZGGRQF, see CGGRQF for further details. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] Q 00070 *> \verbatim 00071 *> Q is COMPLEX*16 array, dimension (LDA,N) 00072 *> The N-by-N unitary matrix Q. 00073 *> \endverbatim 00074 *> 00075 *> \param[out] R 00076 *> \verbatim 00077 *> R is COMPLEX*16 array, dimension (LDA,MAX(M,N)) 00078 *> \endverbatim 00079 *> 00080 *> \param[in] LDA 00081 *> \verbatim 00082 *> LDA is INTEGER 00083 *> The leading dimension of the arrays A, AF, R and Q. 00084 *> LDA >= max(M,N). 00085 *> \endverbatim 00086 *> 00087 *> \param[out] TAUA 00088 *> \verbatim 00089 *> TAUA is COMPLEX*16 array, dimension (min(M,N)) 00090 *> The scalar factors of the elementary reflectors, as returned 00091 *> by DGGQRC. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] B 00095 *> \verbatim 00096 *> B is COMPLEX*16 array, dimension (LDB,N) 00097 *> On entry, the P-by-N matrix A. 00098 *> \endverbatim 00099 *> 00100 *> \param[out] BF 00101 *> \verbatim 00102 *> BF is COMPLEX*16 array, dimension (LDB,N) 00103 *> Details of the GQR factorization of A and B, as returned 00104 *> by ZGGRQF, see CGGRQF for further details. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Z 00108 *> \verbatim 00109 *> Z is DOUBLE PRECISION array, dimension (LDB,P) 00110 *> The P-by-P unitary matrix Z. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] T 00114 *> \verbatim 00115 *> T is COMPLEX*16 array, dimension (LDB,max(P,N)) 00116 *> \endverbatim 00117 *> 00118 *> \param[out] BWK 00119 *> \verbatim 00120 *> BWK is COMPLEX*16 array, dimension (LDB,N) 00121 *> \endverbatim 00122 *> 00123 *> \param[in] LDB 00124 *> \verbatim 00125 *> LDB is INTEGER 00126 *> The leading dimension of the arrays B, BF, Z and T. 00127 *> LDB >= max(P,N). 00128 *> \endverbatim 00129 *> 00130 *> \param[out] TAUB 00131 *> \verbatim 00132 *> TAUB is COMPLEX*16 array, dimension (min(P,N)) 00133 *> The scalar factors of the elementary reflectors, as returned 00134 *> by DGGRQF. 00135 *> \endverbatim 00136 *> 00137 *> \param[out] WORK 00138 *> \verbatim 00139 *> WORK is COMPLEX*16 array, dimension (LWORK) 00140 *> \endverbatim 00141 *> 00142 *> \param[in] LWORK 00143 *> \verbatim 00144 *> LWORK is INTEGER 00145 *> The dimension of the array WORK, LWORK >= max(M,P,N)**2. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] RWORK 00149 *> \verbatim 00150 *> RWORK is DOUBLE PRECISION array, dimension (M) 00151 *> \endverbatim 00152 *> 00153 *> \param[out] RESULT 00154 *> \verbatim 00155 *> RESULT is DOUBLE PRECISION array, dimension (4) 00156 *> The test ratios: 00157 *> RESULT(1) = norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP) 00158 *> RESULT(2) = norm( T*Q - Z'*B ) / (MAX(P,N)*norm(B)*ULP) 00159 *> RESULT(3) = norm( I - Q'*Q ) / ( N*ULP ) 00160 *> RESULT(4) = norm( I - Z'*Z ) / ( P*ULP ) 00161 *> \endverbatim 00162 * 00163 * Authors: 00164 * ======== 00165 * 00166 *> \author Univ. of Tennessee 00167 *> \author Univ. of California Berkeley 00168 *> \author Univ. of Colorado Denver 00169 *> \author NAG Ltd. 00170 * 00171 *> \date November 2011 00172 * 00173 *> \ingroup complex16_eig 00174 * 00175 * ===================================================================== 00176 SUBROUTINE ZGRQTS( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, 00177 $ BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT ) 00178 * 00179 * -- LAPACK test routine (version 3.4.0) -- 00180 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00182 * November 2011 00183 * 00184 * .. Scalar Arguments .. 00185 INTEGER LDA, LDB, LWORK, M, N, P 00186 * .. 00187 * .. Array Arguments .. 00188 DOUBLE PRECISION RESULT( 4 ), RWORK( * ) 00189 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00190 $ BF( LDB, * ), BWK( LDB, * ), Q( LDA, * ), 00191 $ R( LDA, * ), T( LDB, * ), TAUA( * ), TAUB( * ), 00192 $ WORK( LWORK ), Z( LDB, * ) 00193 * .. 00194 * 00195 * ===================================================================== 00196 * 00197 * .. Parameters .. 00198 DOUBLE PRECISION ZERO, ONE 00199 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00200 COMPLEX*16 CZERO, CONE 00201 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00202 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00203 COMPLEX*16 CROGUE 00204 PARAMETER ( CROGUE = ( -1.0D+10, 0.0D+0 ) ) 00205 * .. 00206 * .. Local Scalars .. 00207 INTEGER INFO 00208 DOUBLE PRECISION ANORM, BNORM, RESID, ULP, UNFL 00209 * .. 00210 * .. External Functions .. 00211 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE 00212 EXTERNAL DLAMCH, ZLANGE, ZLANHE 00213 * .. 00214 * .. External Subroutines .. 00215 EXTERNAL ZGEMM, ZGGRQF, ZHERK, ZLACPY, ZLASET, ZUNGQR, 00216 $ ZUNGRQ 00217 * .. 00218 * .. Intrinsic Functions .. 00219 INTRINSIC DBLE, MAX, MIN 00220 * .. 00221 * .. Executable Statements .. 00222 * 00223 ULP = DLAMCH( 'Precision' ) 00224 UNFL = DLAMCH( 'Safe minimum' ) 00225 * 00226 * Copy the matrix A to the array AF. 00227 * 00228 CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) 00229 CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB ) 00230 * 00231 ANORM = MAX( ZLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 00232 BNORM = MAX( ZLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 00233 * 00234 * Factorize the matrices A and B in the arrays AF and BF. 00235 * 00236 CALL ZGGRQF( M, P, N, AF, LDA, TAUA, BF, LDB, TAUB, WORK, LWORK, 00237 $ INFO ) 00238 * 00239 * Generate the N-by-N matrix Q 00240 * 00241 CALL ZLASET( 'Full', N, N, CROGUE, CROGUE, Q, LDA ) 00242 IF( M.LE.N ) THEN 00243 IF( M.GT.0 .AND. M.LT.N ) 00244 $ CALL ZLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA ) 00245 IF( M.GT.1 ) 00246 $ CALL ZLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA, 00247 $ Q( N-M+2, N-M+1 ), LDA ) 00248 ELSE 00249 IF( N.GT.1 ) 00250 $ CALL ZLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA, 00251 $ Q( 2, 1 ), LDA ) 00252 END IF 00253 CALL ZUNGRQ( N, N, MIN( M, N ), Q, LDA, TAUA, WORK, LWORK, INFO ) 00254 * 00255 * Generate the P-by-P matrix Z 00256 * 00257 CALL ZLASET( 'Full', P, P, CROGUE, CROGUE, Z, LDB ) 00258 IF( P.GT.1 ) 00259 $ CALL ZLACPY( 'Lower', P-1, N, BF( 2, 1 ), LDB, Z( 2, 1 ), LDB ) 00260 CALL ZUNGQR( P, P, MIN( P, N ), Z, LDB, TAUB, WORK, LWORK, INFO ) 00261 * 00262 * Copy R 00263 * 00264 CALL ZLASET( 'Full', M, N, CZERO, CZERO, R, LDA ) 00265 IF( M.LE.N ) THEN 00266 CALL ZLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, R( 1, N-M+1 ), 00267 $ LDA ) 00268 ELSE 00269 CALL ZLACPY( 'Full', M-N, N, AF, LDA, R, LDA ) 00270 CALL ZLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, R( M-N+1, 1 ), 00271 $ LDA ) 00272 END IF 00273 * 00274 * Copy T 00275 * 00276 CALL ZLASET( 'Full', P, N, CZERO, CZERO, T, LDB ) 00277 CALL ZLACPY( 'Upper', P, N, BF, LDB, T, LDB ) 00278 * 00279 * Compute R - A*Q' 00280 * 00281 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M, N, N, -CONE, 00282 $ A, LDA, Q, LDA, CONE, R, LDA ) 00283 * 00284 * Compute norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP ) . 00285 * 00286 RESID = ZLANGE( '1', M, N, R, LDA, RWORK ) 00287 IF( ANORM.GT.ZERO ) THEN 00288 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / 00289 $ ULP 00290 ELSE 00291 RESULT( 1 ) = ZERO 00292 END IF 00293 * 00294 * Compute T*Q - Z'*B 00295 * 00296 CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE, 00297 $ Z, LDB, B, LDB, CZERO, BWK, LDB ) 00298 CALL ZGEMM( 'No transpose', 'No transpose', P, N, N, CONE, T, LDB, 00299 $ Q, LDA, -CONE, BWK, LDB ) 00300 * 00301 * Compute norm( T*Q - Z'*B ) / ( MAX(P,N)*norm(A)*ULP ) . 00302 * 00303 RESID = ZLANGE( '1', P, N, BWK, LDB, RWORK ) 00304 IF( BNORM.GT.ZERO ) THEN 00305 RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, M ) ) ) / BNORM ) / 00306 $ ULP 00307 ELSE 00308 RESULT( 2 ) = ZERO 00309 END IF 00310 * 00311 * Compute I - Q*Q' 00312 * 00313 CALL ZLASET( 'Full', N, N, CZERO, CONE, R, LDA ) 00314 CALL ZHERK( 'Upper', 'No Transpose', N, N, -ONE, Q, LDA, ONE, R, 00315 $ LDA ) 00316 * 00317 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00318 * 00319 RESID = ZLANHE( '1', 'Upper', N, R, LDA, RWORK ) 00320 RESULT( 3 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP 00321 * 00322 * Compute I - Z'*Z 00323 * 00324 CALL ZLASET( 'Full', P, P, CZERO, CONE, T, LDB ) 00325 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, Z, LDB, 00326 $ ONE, T, LDB ) 00327 * 00328 * Compute norm( I - Z'*Z ) / ( P*ULP ) . 00329 * 00330 RESID = ZLANHE( '1', 'Upper', P, T, LDB, RWORK ) 00331 RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP 00332 * 00333 RETURN 00334 * 00335 * End of ZGRQTS 00336 * 00337 END