 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 00002 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 00003 $ LWORK, RWORK, RESULT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IWORK( * ) 00014 DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) 00015 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00016 $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), 00017 $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * ZGSVTS tests ZGGSVD, which computes the GSVD of an M-by-N matrix A 00024 * and a P-by-N matrix B: 00025 * U'*A*Q = D1*R and V'*B*Q = D2*R. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * M (input) INTEGER 00031 * The number of rows of the matrix A. M >= 0. 00032 * 00033 * P (input) INTEGER 00034 * The number of rows of the matrix B. P >= 0. 00035 * 00036 * N (input) INTEGER 00037 * The number of columns of the matrices A and B. N >= 0. 00038 * 00039 * A (input) COMPLEX*16 array, dimension (LDA,M) 00040 * The M-by-N matrix A. 00041 * 00042 * AF (output) COMPLEX*16 array, dimension (LDA,N) 00043 * Details of the GSVD of A and B, as returned by ZGGSVD, 00044 * see ZGGSVD for further details. 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the arrays A and AF. 00048 * LDA >= max( 1,M ). 00049 * 00050 * B (input) COMPLEX*16 array, dimension (LDB,P) 00051 * On entry, the P-by-N matrix B. 00052 * 00053 * BF (output) COMPLEX*16 array, dimension (LDB,N) 00054 * Details of the GSVD of A and B, as returned by ZGGSVD, 00055 * see ZGGSVD for further details. 00056 * 00057 * LDB (input) INTEGER 00058 * The leading dimension of the arrays B and BF. 00059 * LDB >= max(1,P). 00060 * 00061 * U (output) COMPLEX*16 array, dimension(LDU,M) 00062 * The M by M unitary matrix U. 00063 * 00064 * LDU (input) INTEGER 00065 * The leading dimension of the array U. LDU >= max(1,M). 00066 * 00067 * V (output) COMPLEX*16 array, dimension(LDV,M) 00068 * The P by P unitary matrix V. 00069 * 00070 * LDV (input) INTEGER 00071 * The leading dimension of the array V. LDV >= max(1,P). 00072 * 00073 * Q (output) COMPLEX*16 array, dimension(LDQ,N) 00074 * The N by N unitary matrix Q. 00075 * 00076 * LDQ (input) INTEGER 00077 * The leading dimension of the array Q. LDQ >= max(1,N). 00078 * 00079 * ALPHA (output) DOUBLE PRECISION array, dimension (N) 00080 * BETA (output) DOUBLE PRECISION array, dimension (N) 00081 * The generalized singular value pairs of A and B, the 00082 * ``diagonal'' matrices D1 and D2 are constructed from 00083 * ALPHA and BETA, see subroutine ZGGSVD for details. 00084 * 00085 * R (output) COMPLEX*16 array, dimension(LDQ,N) 00086 * The upper triangular matrix R. 00087 * 00088 * LDR (input) INTEGER 00089 * The leading dimension of the array R. LDR >= max(1,N). 00090 * 00091 * IWORK (workspace) INTEGER array, dimension (N) 00092 * 00093 * WORK (workspace) COMPLEX*16 array, dimension (LWORK) 00094 * 00095 * LWORK (input) INTEGER 00096 * The dimension of the array WORK, 00097 * LWORK >= max(M,P,N)*max(M,P,N). 00098 * 00099 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(M,P,N)) 00100 * 00101 * RESULT (output) DOUBLE PRECISION array, dimension (5) 00102 * The test ratios: 00103 * RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) 00104 * RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) 00105 * RESULT(3) = norm( I - U'*U ) / ( M*ULP ) 00106 * RESULT(4) = norm( I - V'*V ) / ( P*ULP ) 00107 * RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) 00108 * RESULT(6) = 0 if ALPHA is in decreasing order; 00109 * = ULPINV otherwise. 00110 * 00111 * ===================================================================== 00112 * 00113 * .. Parameters .. 00114 DOUBLE PRECISION ZERO, ONE 00115 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00116 COMPLEX*16 CZERO, CONE 00117 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00118 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00119 * .. 00120 * .. Local Scalars .. 00121 INTEGER I, INFO, J, K, L 00122 DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL 00123 * .. 00124 * .. External Functions .. 00125 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE 00126 EXTERNAL DLAMCH, ZLANGE, ZLANHE 00127 * .. 00128 * .. External Subroutines .. 00129 EXTERNAL DCOPY, ZGEMM, ZGGSVD, ZHERK, ZLACPY, ZLASET 00130 * .. 00131 * .. Intrinsic Functions .. 00132 INTRINSIC DBLE, MAX, MIN 00133 * .. 00134 * .. Executable Statements .. 00135 * 00136 ULP = DLAMCH( 'Precision' ) 00137 ULPINV = ONE / ULP 00138 UNFL = DLAMCH( 'Safe minimum' ) 00139 * 00140 * Copy the matrix A to the array AF. 00141 * 00142 CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) 00143 CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB ) 00144 * 00145 ANORM = MAX( ZLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 00146 BNORM = MAX( ZLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 00147 * 00148 * Factorize the matrices A and B in the arrays AF and BF. 00149 * 00150 CALL ZGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, 00151 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, RWORK, 00152 $ IWORK, INFO ) 00153 * 00154 * Copy R 00155 * 00156 DO 20 I = 1, MIN( K+L, M ) 00157 DO 10 J = I, K + L 00158 R( I, J ) = AF( I, N-K-L+J ) 00159 10 CONTINUE 00160 20 CONTINUE 00161 * 00162 IF( M-K-L.LT.0 ) THEN 00163 DO 40 I = M + 1, K + L 00164 DO 30 J = I, K + L 00165 R( I, J ) = BF( I-K, N-K-L+J ) 00166 30 CONTINUE 00167 40 CONTINUE 00168 END IF 00169 * 00170 * Compute A:= U'*A*Q - D1*R 00171 * 00172 CALL ZGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA, 00173 $ Q, LDQ, CZERO, WORK, LDA ) 00174 * 00175 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE, 00176 $ U, LDU, WORK, LDA, CZERO, A, LDA ) 00177 * 00178 DO 60 I = 1, K 00179 DO 50 J = I, K + L 00180 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) 00181 50 CONTINUE 00182 60 CONTINUE 00183 * 00184 DO 80 I = K + 1, MIN( K+L, M ) 00185 DO 70 J = I, K + L 00186 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) 00187 70 CONTINUE 00188 80 CONTINUE 00189 * 00190 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . 00191 * 00192 RESID = ZLANGE( '1', M, N, A, LDA, RWORK ) 00193 IF( ANORM.GT.ZERO ) THEN 00194 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / 00195 $ ULP 00196 ELSE 00197 RESULT( 1 ) = ZERO 00198 END IF 00199 * 00200 * Compute B := V'*B*Q - D2*R 00201 * 00202 CALL ZGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB, 00203 $ Q, LDQ, CZERO, WORK, LDB ) 00204 * 00205 CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE, 00206 $ V, LDV, WORK, LDB, CZERO, B, LDB ) 00207 * 00208 DO 100 I = 1, L 00209 DO 90 J = I, L 00210 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) 00211 90 CONTINUE 00212 100 CONTINUE 00213 * 00214 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . 00215 * 00216 RESID = ZLANGE( '1', P, N, B, LDB, RWORK ) 00217 IF( BNORM.GT.ZERO ) THEN 00218 RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, N ) ) ) / BNORM ) / 00219 $ ULP 00220 ELSE 00221 RESULT( 2 ) = ZERO 00222 END IF 00223 * 00224 * Compute I - U'*U 00225 * 00226 CALL ZLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ ) 00227 CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU, 00228 $ ONE, WORK, LDU ) 00229 * 00230 * Compute norm( I - U'*U ) / ( M * ULP ) . 00231 * 00232 RESID = ZLANHE( '1', 'Upper', M, WORK, LDU, RWORK ) 00233 RESULT( 3 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / ULP 00234 * 00235 * Compute I - V'*V 00236 * 00237 CALL ZLASET( 'Full', P, P, CZERO, CONE, WORK, LDV ) 00238 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV, 00239 $ ONE, WORK, LDV ) 00240 * 00241 * Compute norm( I - V'*V ) / ( P * ULP ) . 00242 * 00243 RESID = ZLANHE( '1', 'Upper', P, WORK, LDV, RWORK ) 00244 RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP 00245 * 00246 * Compute I - Q'*Q 00247 * 00248 CALL ZLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ ) 00249 CALL ZHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ, 00250 $ ONE, WORK, LDQ ) 00251 * 00252 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00253 * 00254 RESID = ZLANHE( '1', 'Upper', N, WORK, LDQ, RWORK ) 00255 RESULT( 5 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP 00256 * 00257 * Check sorting 00258 * 00259 CALL DCOPY( N, ALPHA, 1, RWORK, 1 ) 00260 DO 110 I = K + 1, MIN( K+L, M ) 00261 J = IWORK( I ) 00262 IF( I.NE.J ) THEN 00263 TEMP = RWORK( I ) 00264 RWORK( I ) = RWORK( J ) 00265 RWORK( J ) = TEMP 00266 END IF 00267 110 CONTINUE 00268 * 00269 RESULT( 6 ) = ZERO 00270 DO 120 I = K + 1, MIN( K+L, M ) - 1 00271 IF( RWORK( I ).LT.RWORK( I+1 ) ) 00272 $ RESULT( 6 ) = ULPINV 00273 120 CONTINUE 00274 * 00275 RETURN 00276 * 00277 * End of ZGSVTS 00278 * 00279 END
1.7.3 
