![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGQRTS 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 ZGQRTS( N, M, P, 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 *> ZGQRTS tests ZGGQRF, which computes the GQR factorization of an 00032 *> N-by-M matrix A and a N-by-P matrix B: A = Q*R and B = Q*T*Z. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] N 00039 *> \verbatim 00040 *> N is INTEGER 00041 *> The number of rows of the matrices A and B. N >= 0. 00042 *> \endverbatim 00043 *> 00044 *> \param[in] M 00045 *> \verbatim 00046 *> M is INTEGER 00047 *> The number of columns of the matrix A. M >= 0. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] P 00051 *> \verbatim 00052 *> P is INTEGER 00053 *> The number of columns of the matrix B. P >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] A 00057 *> \verbatim 00058 *> A is COMPLEX*16 array, dimension (LDA,M) 00059 *> The N-by-M matrix A. 00060 *> \endverbatim 00061 *> 00062 *> \param[out] AF 00063 *> \verbatim 00064 *> AF is COMPLEX*16 array, dimension (LDA,N) 00065 *> Details of the GQR factorization of A and B, as returned 00066 *> by ZGGQRF, see CGGQRF for further details. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] Q 00070 *> \verbatim 00071 *> Q is COMPLEX*16 array, dimension (LDA,N) 00072 *> The M-by-M 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 ZGGQRF. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] B 00095 *> \verbatim 00096 *> B is COMPLEX*16 array, dimension (LDB,P) 00097 *> On entry, the N-by-P 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 ZGGQRF, see CGGQRF for further details. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Z 00108 *> \verbatim 00109 *> Z is COMPLEX*16 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(N,M,P)**2. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] RWORK 00149 *> \verbatim 00150 *> RWORK is DOUBLE PRECISION array, dimension (max(N,M,P)) 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 - Q'*A ) / ( MAX(M,N)*norm(A)*ULP) 00158 *> RESULT(2) = norm( T*Z - Q'*B ) / (MAX(P,N)*norm(B)*ULP) 00159 *> RESULT(3) = norm( I - Q'*Q ) / ( M*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 ZGQRTS( N, M, P, 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, ZGGQRF, 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', N, M, A, LDA, AF, LDA ) 00229 CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB ) 00230 * 00231 ANORM = MAX( ZLANGE( '1', N, M, A, LDA, RWORK ), UNFL ) 00232 BNORM = MAX( ZLANGE( '1', N, P, B, LDB, RWORK ), UNFL ) 00233 * 00234 * Factorize the matrices A and B in the arrays AF and BF. 00235 * 00236 CALL ZGGQRF( N, M, P, 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 CALL ZLACPY( 'Lower', N-1, M, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) 00243 CALL ZUNGQR( N, N, MIN( N, M ), Q, LDA, TAUA, WORK, LWORK, INFO ) 00244 * 00245 * Generate the P-by-P matrix Z 00246 * 00247 CALL ZLASET( 'Full', P, P, CROGUE, CROGUE, Z, LDB ) 00248 IF( N.LE.P ) THEN 00249 IF( N.GT.0 .AND. N.LT.P ) 00250 $ CALL ZLACPY( 'Full', N, P-N, BF, LDB, Z( P-N+1, 1 ), LDB ) 00251 IF( N.GT.1 ) 00252 $ CALL ZLACPY( 'Lower', N-1, N-1, BF( 2, P-N+1 ), LDB, 00253 $ Z( P-N+2, P-N+1 ), LDB ) 00254 ELSE 00255 IF( P.GT.1 ) 00256 $ CALL ZLACPY( 'Lower', P-1, P-1, BF( N-P+2, 1 ), LDB, 00257 $ Z( 2, 1 ), LDB ) 00258 END IF 00259 CALL ZUNGRQ( P, P, MIN( N, P ), Z, LDB, TAUB, WORK, LWORK, INFO ) 00260 * 00261 * Copy R 00262 * 00263 CALL ZLASET( 'Full', N, M, CZERO, CZERO, R, LDA ) 00264 CALL ZLACPY( 'Upper', N, M, AF, LDA, R, LDA ) 00265 * 00266 * Copy T 00267 * 00268 CALL ZLASET( 'Full', N, P, CZERO, CZERO, T, LDB ) 00269 IF( N.LE.P ) THEN 00270 CALL ZLACPY( 'Upper', N, N, BF( 1, P-N+1 ), LDB, T( 1, P-N+1 ), 00271 $ LDB ) 00272 ELSE 00273 CALL ZLACPY( 'Full', N-P, P, BF, LDB, T, LDB ) 00274 CALL ZLACPY( 'Upper', P, P, BF( N-P+1, 1 ), LDB, T( N-P+1, 1 ), 00275 $ LDB ) 00276 END IF 00277 * 00278 * Compute R - Q'*A 00279 * 00280 CALL ZGEMM( 'Conjugate transpose', 'No transpose', N, M, N, -CONE, 00281 $ Q, LDA, A, LDA, CONE, R, LDA ) 00282 * 00283 * Compute norm( R - Q'*A ) / ( MAX(M,N)*norm(A)*ULP ) . 00284 * 00285 RESID = ZLANGE( '1', N, M, R, LDA, RWORK ) 00286 IF( ANORM.GT.ZERO ) THEN 00287 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / 00288 $ ULP 00289 ELSE 00290 RESULT( 1 ) = ZERO 00291 END IF 00292 * 00293 * Compute T*Z - Q'*B 00294 * 00295 CALL ZGEMM( 'No Transpose', 'No transpose', N, P, P, CONE, T, LDB, 00296 $ Z, LDB, CZERO, BWK, LDB ) 00297 CALL ZGEMM( 'Conjugate transpose', 'No transpose', N, P, N, -CONE, 00298 $ Q, LDA, B, LDB, CONE, BWK, LDB ) 00299 * 00300 * Compute norm( T*Z - Q'*B ) / ( MAX(P,N)*norm(A)*ULP ) . 00301 * 00302 RESID = ZLANGE( '1', N, P, BWK, LDB, RWORK ) 00303 IF( BNORM.GT.ZERO ) THEN 00304 RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, N ) ) ) / BNORM ) / 00305 $ ULP 00306 ELSE 00307 RESULT( 2 ) = ZERO 00308 END IF 00309 * 00310 * Compute I - Q'*Q 00311 * 00312 CALL ZLASET( 'Full', N, N, CZERO, CONE, R, LDA ) 00313 CALL ZHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDA, 00314 $ ONE, R, LDA ) 00315 * 00316 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00317 * 00318 RESID = ZLANHE( '1', 'Upper', N, R, LDA, RWORK ) 00319 RESULT( 3 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP 00320 * 00321 * Compute I - Z'*Z 00322 * 00323 CALL ZLASET( 'Full', P, P, CZERO, CONE, T, LDB ) 00324 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, Z, LDB, 00325 $ ONE, T, LDB ) 00326 * 00327 * Compute norm( I - Z'*Z ) / ( P*ULP ) . 00328 * 00329 RESID = ZLANHE( '1', 'Upper', P, T, LDB, RWORK ) 00330 RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP 00331 * 00332 RETURN 00333 * 00334 * End of ZGQRTS 00335 * 00336 END