![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZGSVTS 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 ZGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 00012 * LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 00013 * LWORK, RWORK, RESULT ) 00014 * 00015 * .. Scalar Arguments .. 00016 * INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IWORK( * ) 00020 * DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) 00021 * COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00022 * $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), 00023 * $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) 00024 * .. 00025 * 00026 * 00027 *> \par Purpose: 00028 * ============= 00029 *> 00030 *> \verbatim 00031 *> 00032 *> ZGSVTS tests ZGGSVD, which computes the GSVD of an M-by-N matrix A 00033 *> and a P-by-N matrix B: 00034 *> U'*A*Q = D1*R and V'*B*Q = D2*R. 00035 *> \endverbatim 00036 * 00037 * Arguments: 00038 * ========== 00039 * 00040 *> \param[in] M 00041 *> \verbatim 00042 *> M is INTEGER 00043 *> The number of rows of the matrix A. M >= 0. 00044 *> \endverbatim 00045 *> 00046 *> \param[in] P 00047 *> \verbatim 00048 *> P is INTEGER 00049 *> The number of rows of the matrix B. P >= 0. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] N 00053 *> \verbatim 00054 *> N is INTEGER 00055 *> The number of columns of the matrices A and B. N >= 0. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] A 00059 *> \verbatim 00060 *> A is COMPLEX*16 array, dimension (LDA,M) 00061 *> The M-by-N matrix A. 00062 *> \endverbatim 00063 *> 00064 *> \param[out] AF 00065 *> \verbatim 00066 *> AF is COMPLEX*16 array, dimension (LDA,N) 00067 *> Details of the GSVD of A and B, as returned by ZGGSVD, 00068 *> see ZGGSVD for further details. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] LDA 00072 *> \verbatim 00073 *> LDA is INTEGER 00074 *> The leading dimension of the arrays A and AF. 00075 *> LDA >= max( 1,M ). 00076 *> \endverbatim 00077 *> 00078 *> \param[in] B 00079 *> \verbatim 00080 *> B is COMPLEX*16 array, dimension (LDB,P) 00081 *> On entry, the P-by-N matrix B. 00082 *> \endverbatim 00083 *> 00084 *> \param[out] BF 00085 *> \verbatim 00086 *> BF is COMPLEX*16 array, dimension (LDB,N) 00087 *> Details of the GSVD of A and B, as returned by ZGGSVD, 00088 *> see ZGGSVD for further details. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LDB 00092 *> \verbatim 00093 *> LDB is INTEGER 00094 *> The leading dimension of the arrays B and BF. 00095 *> LDB >= max(1,P). 00096 *> \endverbatim 00097 *> 00098 *> \param[out] U 00099 *> \verbatim 00100 *> U is COMPLEX*16 array, dimension(LDU,M) 00101 *> The M by M unitary matrix U. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] LDU 00105 *> \verbatim 00106 *> LDU is INTEGER 00107 *> The leading dimension of the array U. LDU >= max(1,M). 00108 *> \endverbatim 00109 *> 00110 *> \param[out] V 00111 *> \verbatim 00112 *> V is COMPLEX*16 array, dimension(LDV,M) 00113 *> The P by P unitary matrix V. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDV 00117 *> \verbatim 00118 *> LDV is INTEGER 00119 *> The leading dimension of the array V. LDV >= max(1,P). 00120 *> \endverbatim 00121 *> 00122 *> \param[out] Q 00123 *> \verbatim 00124 *> Q is COMPLEX*16 array, dimension(LDQ,N) 00125 *> The N by N unitary matrix Q. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] LDQ 00129 *> \verbatim 00130 *> LDQ is INTEGER 00131 *> The leading dimension of the array Q. LDQ >= max(1,N). 00132 *> \endverbatim 00133 *> 00134 *> \param[out] ALPHA 00135 *> \verbatim 00136 *> ALPHA is DOUBLE PRECISION array, dimension (N) 00137 *> \endverbatim 00138 *> 00139 *> \param[out] BETA 00140 *> \verbatim 00141 *> BETA is DOUBLE PRECISION array, dimension (N) 00142 *> 00143 *> The generalized singular value pairs of A and B, the 00144 *> ``diagonal'' matrices D1 and D2 are constructed from 00145 *> ALPHA and BETA, see subroutine ZGGSVD for details. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] R 00149 *> \verbatim 00150 *> R is COMPLEX*16 array, dimension(LDQ,N) 00151 *> The upper triangular matrix R. 00152 *> \endverbatim 00153 *> 00154 *> \param[in] LDR 00155 *> \verbatim 00156 *> LDR is INTEGER 00157 *> The leading dimension of the array R. LDR >= max(1,N). 00158 *> \endverbatim 00159 *> 00160 *> \param[out] IWORK 00161 *> \verbatim 00162 *> IWORK is INTEGER array, dimension (N) 00163 *> \endverbatim 00164 *> 00165 *> \param[out] WORK 00166 *> \verbatim 00167 *> WORK is COMPLEX*16 array, dimension (LWORK) 00168 *> \endverbatim 00169 *> 00170 *> \param[in] LWORK 00171 *> \verbatim 00172 *> LWORK is INTEGER 00173 *> The dimension of the array WORK, 00174 *> LWORK >= max(M,P,N)*max(M,P,N). 00175 *> \endverbatim 00176 *> 00177 *> \param[out] RWORK 00178 *> \verbatim 00179 *> RWORK is DOUBLE PRECISION array, dimension (max(M,P,N)) 00180 *> \endverbatim 00181 *> 00182 *> \param[out] RESULT 00183 *> \verbatim 00184 *> RESULT is DOUBLE PRECISION array, dimension (5) 00185 *> The test ratios: 00186 *> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) 00187 *> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) 00188 *> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) 00189 *> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) 00190 *> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) 00191 *> RESULT(6) = 0 if ALPHA is in decreasing order; 00192 *> = ULPINV otherwise. 00193 *> \endverbatim 00194 * 00195 * Authors: 00196 * ======== 00197 * 00198 *> \author Univ. of Tennessee 00199 *> \author Univ. of California Berkeley 00200 *> \author Univ. of Colorado Denver 00201 *> \author NAG Ltd. 00202 * 00203 *> \date November 2011 00204 * 00205 *> \ingroup complex16_eig 00206 * 00207 * ===================================================================== 00208 SUBROUTINE ZGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 00209 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 00210 $ LWORK, RWORK, RESULT ) 00211 * 00212 * -- LAPACK test routine (version 3.4.0) -- 00213 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00215 * November 2011 00216 * 00217 * .. Scalar Arguments .. 00218 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 00219 * .. 00220 * .. Array Arguments .. 00221 INTEGER IWORK( * ) 00222 DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) 00223 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00224 $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), 00225 $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) 00226 * .. 00227 * 00228 * ===================================================================== 00229 * 00230 * .. Parameters .. 00231 DOUBLE PRECISION ZERO, ONE 00232 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00233 COMPLEX*16 CZERO, CONE 00234 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00235 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00236 * .. 00237 * .. Local Scalars .. 00238 INTEGER I, INFO, J, K, L 00239 DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL 00240 * .. 00241 * .. External Functions .. 00242 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE 00243 EXTERNAL DLAMCH, ZLANGE, ZLANHE 00244 * .. 00245 * .. External Subroutines .. 00246 EXTERNAL DCOPY, ZGEMM, ZGGSVD, ZHERK, ZLACPY, ZLASET 00247 * .. 00248 * .. Intrinsic Functions .. 00249 INTRINSIC DBLE, MAX, MIN 00250 * .. 00251 * .. Executable Statements .. 00252 * 00253 ULP = DLAMCH( 'Precision' ) 00254 ULPINV = ONE / ULP 00255 UNFL = DLAMCH( 'Safe minimum' ) 00256 * 00257 * Copy the matrix A to the array AF. 00258 * 00259 CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) 00260 CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB ) 00261 * 00262 ANORM = MAX( ZLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 00263 BNORM = MAX( ZLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 00264 * 00265 * Factorize the matrices A and B in the arrays AF and BF. 00266 * 00267 CALL ZGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, 00268 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, 00269 $ IWORK, INFO ) 00270 * 00271 * Copy R 00272 * 00273 DO 20 I = 1, MIN( K+L, M ) 00274 DO 10 J = I, K + L 00275 R( I, J ) = AF( I, N-K-L+J ) 00276 10 CONTINUE 00277 20 CONTINUE 00278 * 00279 IF( M-K-L.LT.0 ) THEN 00280 DO 40 I = M + 1, K + L 00281 DO 30 J = I, K + L 00282 R( I, J ) = BF( I-K, N-K-L+J ) 00283 30 CONTINUE 00284 40 CONTINUE 00285 END IF 00286 * 00287 * Compute A:= U'*A*Q - D1*R 00288 * 00289 CALL ZGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA, 00290 $ Q, LDQ, CZERO, WORK, LDA ) 00291 * 00292 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE, 00293 $ U, LDU, WORK, LDA, CZERO, A, LDA ) 00294 * 00295 DO 60 I = 1, K 00296 DO 50 J = I, K + L 00297 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) 00298 50 CONTINUE 00299 60 CONTINUE 00300 * 00301 DO 80 I = K + 1, MIN( K+L, M ) 00302 DO 70 J = I, K + L 00303 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) 00304 70 CONTINUE 00305 80 CONTINUE 00306 * 00307 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . 00308 * 00309 RESID = ZLANGE( '1', M, N, A, LDA, RWORK ) 00310 IF( ANORM.GT.ZERO ) THEN 00311 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / 00312 $ ULP 00313 ELSE 00314 RESULT( 1 ) = ZERO 00315 END IF 00316 * 00317 * Compute B := V'*B*Q - D2*R 00318 * 00319 CALL ZGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB, 00320 $ Q, LDQ, CZERO, WORK, LDB ) 00321 * 00322 CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE, 00323 $ V, LDV, WORK, LDB, CZERO, B, LDB ) 00324 * 00325 DO 100 I = 1, L 00326 DO 90 J = I, L 00327 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) 00328 90 CONTINUE 00329 100 CONTINUE 00330 * 00331 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . 00332 * 00333 RESID = ZLANGE( '1', P, N, B, LDB, RWORK ) 00334 IF( BNORM.GT.ZERO ) THEN 00335 RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, N ) ) ) / BNORM ) / 00336 $ ULP 00337 ELSE 00338 RESULT( 2 ) = ZERO 00339 END IF 00340 * 00341 * Compute I - U'*U 00342 * 00343 CALL ZLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ ) 00344 CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU, 00345 $ ONE, WORK, LDU ) 00346 * 00347 * Compute norm( I - U'*U ) / ( M * ULP ) . 00348 * 00349 RESID = ZLANHE( '1', 'Upper', M, WORK, LDU, RWORK ) 00350 RESULT( 3 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / ULP 00351 * 00352 * Compute I - V'*V 00353 * 00354 CALL ZLASET( 'Full', P, P, CZERO, CONE, WORK, LDV ) 00355 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV, 00356 $ ONE, WORK, LDV ) 00357 * 00358 * Compute norm( I - V'*V ) / ( P * ULP ) . 00359 * 00360 RESID = ZLANHE( '1', 'Upper', P, WORK, LDV, RWORK ) 00361 RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP 00362 * 00363 * Compute I - Q'*Q 00364 * 00365 CALL ZLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ ) 00366 CALL ZHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ, 00367 $ ONE, WORK, LDQ ) 00368 * 00369 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00370 * 00371 RESID = ZLANHE( '1', 'Upper', N, WORK, LDQ, RWORK ) 00372 RESULT( 5 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP 00373 * 00374 * Check sorting 00375 * 00376 CALL DCOPY( N, ALPHA, 1, RWORK, 1 ) 00377 DO 110 I = K + 1, MIN( K+L, M ) 00378 J = IWORK( I ) 00379 IF( I.NE.J ) THEN 00380 TEMP = RWORK( I ) 00381 RWORK( I ) = RWORK( J ) 00382 RWORK( J ) = TEMP 00383 END IF 00384 110 CONTINUE 00385 * 00386 RESULT( 6 ) = ZERO 00387 DO 120 I = K + 1, MIN( K+L, M ) - 1 00388 IF( RWORK( I ).LT.RWORK( I+1 ) ) 00389 $ RESULT( 6 ) = ULPINV 00390 120 CONTINUE 00391 * 00392 RETURN 00393 * 00394 * End of ZGSVTS 00395 * 00396 END