|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, 00002 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, 00003 $ RWORK, RESULT ) 00004 IMPLICIT NONE 00005 * 00006 * Originally xGSVTS 00007 * -- LAPACK test routine (version 3.3.0) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * November 2010 00010 * 00011 * Adapted to CCSDTS by 00012 * July 2010 00013 * 00014 * .. Scalar Arguments .. 00015 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q 00016 * .. 00017 * .. Array Arguments .. 00018 INTEGER IWORK( * ) 00019 REAL RESULT( 9 ), RWORK( * ), THETA( * ) 00020 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00021 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ), 00022 $ XF( LDX, * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * CCSDTS tests CUNCSD, which, given an M-by-M partitioned unitary 00029 * matrix X, 00030 * Q M-Q 00031 * X = [ X11 X12 ] P , 00032 * [ X21 X22 ] M-P 00033 * 00034 * computes the CSD 00035 * 00036 * [ U1 ]**T * [ X11 X12 ] * [ V1 ] 00037 * [ U2 ] [ X21 X22 ] [ V2 ] 00038 * 00039 * [ I 0 0 | 0 0 0 ] 00040 * [ 0 C 0 | 0 -S 0 ] 00041 * [ 0 0 0 | 0 0 -I ] 00042 * = [---------------------] = [ D11 D12 ] . 00043 * [ 0 0 0 | I 0 0 ] [ D21 D22 ] 00044 * [ 0 S 0 | 0 C 0 ] 00045 * [ 0 0 I | 0 0 0 ] 00046 * 00047 * Arguments 00048 * ========= 00049 * 00050 * M (input) INTEGER 00051 * The number of rows of the matrix X. M >= 0. 00052 * 00053 * P (input) INTEGER 00054 * The number of rows of the matrix X11. P >= 0. 00055 * 00056 * Q (input) INTEGER 00057 * The number of columns of the matrix X11. Q >= 0. 00058 * 00059 * X (input) COMPLEX array, dimension (LDX,M) 00060 * The M-by-M matrix X. 00061 * 00062 * XF (output) COMPLEX array, dimension (LDX,M) 00063 * Details of the CSD of X, as returned by CUNCSD; 00064 * see CUNCSD for further details. 00065 * 00066 * LDX (input) INTEGER 00067 * The leading dimension of the arrays X and XF. 00068 * LDX >= max( 1,M ). 00069 * 00070 * U1 (output) COMPLEX array, dimension(LDU1,P) 00071 * The P-by-P unitary matrix U1. 00072 * 00073 * LDU1 (input) INTEGER 00074 * The leading dimension of the array U1. LDU >= max(1,P). 00075 * 00076 * U2 (output) COMPLEX array, dimension(LDU2,M-P) 00077 * The (M-P)-by-(M-P) unitary matrix U2. 00078 * 00079 * LDU2 (input) INTEGER 00080 * The leading dimension of the array U2. LDU >= max(1,M-P). 00081 * 00082 * V1T (output) COMPLEX array, dimension(LDV1T,Q) 00083 * The Q-by-Q unitary matrix V1T. 00084 * 00085 * LDV1T (input) INTEGER 00086 * The leading dimension of the array V1T. LDV1T >= 00087 * max(1,Q). 00088 * 00089 * V2T (output) COMPLEX array, dimension(LDV2T,M-Q) 00090 * The (M-Q)-by-(M-Q) unitary matrix V2T. 00091 * 00092 * LDV2T (input) INTEGER 00093 * The leading dimension of the array V2T. LDV2T >= 00094 * max(1,M-Q). 00095 * 00096 * THETA (output) REAL array, dimension MIN(P,M-P,Q,M-Q) 00097 * The CS values of X; the essentially diagonal matrices C and 00098 * S are constructed from THETA; see subroutine CUNCSD for 00099 * details. 00100 * 00101 * IWORK (workspace) INTEGER array, dimension (M) 00102 * 00103 * WORK (workspace) COMPLEX array, dimension (LWORK) 00104 * 00105 * LWORK (input) INTEGER 00106 * The dimension of the array WORK 00107 * 00108 * RWORK (workspace) REAL array 00109 * 00110 * RESULT (output) REAL array, dimension (9) 00111 * The test ratios: 00112 * RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) 00113 * RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) 00114 * RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) 00115 * RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) 00116 * RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP ) 00117 * RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP ) 00118 * RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP ) 00119 * RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP ) 00120 * RESULT(9) = 0 if THETA is in increasing order and 00121 * all angles are in [0,pi/2]; 00122 * = ULPINV otherwise. 00123 * ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). ) 00124 * 00125 * ===================================================================== 00126 * 00127 * .. Parameters .. 00128 REAL PIOVER2, REALONE, REALZERO 00129 PARAMETER ( PIOVER2 = 1.57079632679489662E0, 00130 $ REALONE = 1.0E0, REALZERO = 0.0E0 ) 00131 COMPLEX ZERO, ONE 00132 PARAMETER ( ZERO = (0.0E0,0.0E0), ONE = (1.0E0,0.0E0) ) 00133 * .. 00134 * .. Local Scalars .. 00135 INTEGER I, INFO, R 00136 REAL EPS2, RESID, ULP, ULPINV 00137 * .. 00138 * .. External Functions .. 00139 REAL SLAMCH, CLANGE, CLANHE 00140 EXTERNAL SLAMCH, CLANGE, CLANHE 00141 * .. 00142 * .. External Subroutines .. 00143 EXTERNAL CGEMM, CLACPY, CLASET, CUNCSD, CHERK 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC REAL, MAX, MIN 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 ULP = SLAMCH( 'Precision' ) 00151 ULPINV = REALONE / ULP 00152 CALL CLASET( 'Full', M, M, ZERO, ONE, WORK, LDX ) 00153 CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX, 00154 $ ONE, WORK, LDX ) 00155 EPS2 = MAX( ULP, 00156 $ CLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) ) 00157 R = MIN( P, M-P, Q, M-Q ) 00158 * 00159 * Copy the matrix X to the array XF. 00160 * 00161 CALL CLACPY( 'Full', M, M, X, LDX, XF, LDX ) 00162 * 00163 * Compute the CSD 00164 * 00165 CALL CUNCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX, 00166 $ XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX, 00167 $ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 00168 $ WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO ) 00169 * 00170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22] 00171 * 00172 CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE, 00173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00174 * 00175 CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE, 00176 $ U1, LDU1, WORK, LDX, ZERO, X, LDX ) 00177 * 00178 DO I = 1, MIN(P,Q)-R 00179 X(I,I) = X(I,I) - ONE 00180 END DO 00181 DO I = 1, R 00182 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = 00183 $ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - CMPLX( COS(THETA(I)), 00184 $ 0.0E0 ) 00185 END DO 00186 * 00187 CALL CGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q, 00188 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00189 * 00190 CALL CGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P, 00191 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX ) 00192 * 00193 DO I = 1, MIN(P,M-Q)-R 00194 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE 00195 END DO 00196 DO I = 1, R 00197 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) = 00198 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) + 00199 $ CMPLX( SIN(THETA(R-I+1)), 0.0E0 ) 00200 END DO 00201 * 00202 CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE, 00203 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00204 * 00205 CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P, 00206 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX ) 00207 * 00208 DO I = 1, MIN(M-P,Q)-R 00209 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE 00210 END DO 00211 DO I = 1, R 00212 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) = 00213 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) - 00214 $ CMPLX( SIN(THETA(R-I+1)), 0.0E0 ) 00215 END DO 00216 * 00217 CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q, 00218 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00219 * 00220 CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P, 00221 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX ) 00222 * 00223 DO I = 1, MIN(M-P,M-Q)-R 00224 X(P+I,Q+I) = X(P+I,Q+I) - ONE 00225 END DO 00226 DO I = 1, R 00227 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) = 00228 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) - 00229 $ CMPLX( COS(THETA(I)), 0.0E0 ) 00230 END DO 00231 * 00232 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) . 00233 * 00234 RESID = CLANGE( '1', P, Q, X, LDX, RWORK ) 00235 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2 00236 * 00237 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) . 00238 * 00239 RESID = CLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK ) 00240 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2 00241 * 00242 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) . 00243 * 00244 RESID = CLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK ) 00245 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2 00246 * 00247 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) . 00248 * 00249 RESID = CLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK ) 00250 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2 00251 * 00252 * Compute I - U1'*U1 00253 * 00254 CALL CLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 ) 00255 CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1, 00256 $ ONE, WORK, LDU1 ) 00257 * 00258 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) . 00259 * 00260 RESID = CLANHE( '1', 'Upper', P, WORK, LDU1, RWORK ) 00261 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP 00262 * 00263 * Compute I - U2'*U2 00264 * 00265 CALL CLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 ) 00266 CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2, 00267 $ LDU2, ONE, WORK, LDU2 ) 00268 * 00269 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) . 00270 * 00271 RESID = CLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK ) 00272 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP 00273 * 00274 * Compute I - V1T*V1T' 00275 * 00276 CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T ) 00277 CALL CHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE, 00278 $ WORK, LDV1T ) 00279 * 00280 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) . 00281 * 00282 RESID = CLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK ) 00283 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP 00284 * 00285 * Compute I - V2T*V2T' 00286 * 00287 CALL CLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T ) 00288 CALL CHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T, 00289 $ ONE, WORK, LDV2T ) 00290 * 00291 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) . 00292 * 00293 RESID = CLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK ) 00294 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP 00295 * 00296 * Check sorting 00297 * 00298 RESULT(9) = REALZERO 00299 DO I = 1, R 00300 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN 00301 RESULT(9) = ULPINV 00302 END IF 00303 IF( I.GT.1) THEN 00304 IF ( THETA(I).LT.THETA(I-1) ) THEN 00305 RESULT(9) = ULPINV 00306 END IF 00307 END IF 00308 END DO 00309 * 00310 RETURN 00311 * 00312 * End of CCSDTS 00313 * 00314 END 00315
1.7.3 
