|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DCSDTS( 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.1) -- 00008 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00009 * November 2006 00010 * 00011 * Adapted to DCSDTS 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 DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * ) 00020 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00021 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ), 00022 $ XF( LDX, * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal 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) DOUBLE PRECISION array, dimension (LDX,M) 00060 * The M-by-M matrix X. 00061 * 00062 * XF (output) DOUBLE PRECISION array, dimension (LDX,M) 00063 * Details of the CSD of X, as returned by DORCSD; 00064 * see DORCSD 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) DOUBLE PRECISION array, dimension(LDU1,P) 00071 * The P-by-P orthogonal matrix U1. 00072 * 00073 * LDU1 (input) INTEGER 00074 * The leading dimension of the array U1. LDU >= max(1,P). 00075 * 00076 * U2 (output) DOUBLE PRECISION array, dimension(LDU2,M-P) 00077 * The (M-P)-by-(M-P) orthogonal matrix U2. 00078 * 00079 * LDU2 (input) INTEGER 00080 * The leading dimension of the array U2. LDU >= max(1,M-P). 00081 * 00082 * V1T (output) DOUBLE PRECISION array, dimension(LDV1T,Q) 00083 * The Q-by-Q orthogonal matrix V1T. 00084 * 00085 * LDV1T (input) INTEGER 00086 * The leading dimension of the array V1T. LDV1T >= 00087 * max(1,Q). 00088 * 00089 * V2T (output) DOUBLE PRECISION array, dimension(LDV2T,M-Q) 00090 * The (M-Q)-by-(M-Q) orthogonal 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) DOUBLE PRECISION 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 DORCSD for 00099 * details. 00100 * 00101 * IWORK (workspace) INTEGER array, dimension (M) 00102 * 00103 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK) 00104 * 00105 * LWORK (input) INTEGER 00106 * The dimension of the array WORK 00107 * 00108 * RWORK (workspace) DOUBLE PRECISION array 00109 * 00110 * RESULT (output) DOUBLE PRECISION 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 DOUBLE PRECISION PIOVER2, REALONE, REALZERO 00129 PARAMETER ( PIOVER2 = 1.57079632679489662D0, 00130 $ REALONE = 1.0D0, REALZERO = 0.0D0 ) 00131 DOUBLE PRECISION ZERO, ONE 00132 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00133 * .. 00134 * .. Local Scalars .. 00135 INTEGER I, INFO, R 00136 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV 00137 * .. 00138 * .. External Functions .. 00139 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00140 EXTERNAL DLAMCH, DLANGE, DLANSY 00141 * .. 00142 * .. External Subroutines .. 00143 EXTERNAL DGEMM, DLACPY, DLASET, DORCSD, DSYRK 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC DBLE, MAX, MIN 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 ULP = DLAMCH( 'Precision' ) 00151 ULPINV = REALONE / ULP 00152 CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDX ) 00153 CALL DSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX, 00154 $ ONE, WORK, LDX ) 00155 EPS2 = MAX( ULP, 00156 $ DLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) ) 00157 R = MIN( P, M-P, Q, M-Q ) 00158 * 00159 * Copy the matrix X to the array XF. 00160 * 00161 CALL DLACPY( 'Full', M, M, X, LDX, XF, LDX ) 00162 * 00163 * Compute the CSD 00164 * 00165 CALL DORCSD( '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, IWORK, INFO ) 00169 * 00170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22] 00171 * 00172 CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE, 00173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00174 * 00175 CALL DGEMM( '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) - COS(THETA(I)) 00184 END DO 00185 * 00186 CALL DGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q, 00187 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00188 * 00189 CALL DGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P, 00190 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX ) 00191 * 00192 DO I = 1, MIN(P,M-Q)-R 00193 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE 00194 END DO 00195 DO I = 1, R 00196 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) = 00197 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) + 00198 $ SIN(THETA(R-I+1)) 00199 END DO 00200 * 00201 CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE, 00202 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00203 * 00204 CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P, 00205 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX ) 00206 * 00207 DO I = 1, MIN(M-P,Q)-R 00208 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE 00209 END DO 00210 DO I = 1, R 00211 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) = 00212 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) - 00213 $ SIN(THETA(R-I+1)) 00214 END DO 00215 * 00216 CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q, 00217 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00218 * 00219 CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P, 00220 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX ) 00221 * 00222 DO I = 1, MIN(M-P,M-Q)-R 00223 X(P+I,Q+I) = X(P+I,Q+I) - ONE 00224 END DO 00225 DO I = 1, R 00226 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) = 00227 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) - 00228 $ COS(THETA(I)) 00229 END DO 00230 * 00231 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) . 00232 * 00233 RESID = DLANGE( '1', P, Q, X, LDX, RWORK ) 00234 RESULT( 1 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2 00235 * 00236 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) . 00237 * 00238 RESID = DLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK ) 00239 RESULT( 2 ) = ( RESID / DBLE(MAX(1,P,M-Q)) ) / EPS2 00240 * 00241 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) . 00242 * 00243 RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK ) 00244 RESULT( 3 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2 00245 * 00246 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) . 00247 * 00248 RESID = DLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK ) 00249 RESULT( 4 ) = ( RESID / DBLE(MAX(1,M-P,M-Q)) ) / EPS2 00250 * 00251 * Compute I - U1'*U1 00252 * 00253 CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 ) 00254 CALL DSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1, 00255 $ ONE, WORK, LDU1 ) 00256 * 00257 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) . 00258 * 00259 RESID = DLANSY( '1', 'Upper', P, WORK, LDU1, RWORK ) 00260 RESULT( 5 ) = ( RESID / DBLE(MAX(1,P)) ) / ULP 00261 * 00262 * Compute I - U2'*U2 00263 * 00264 CALL DLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 ) 00265 CALL DSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2, 00266 $ LDU2, ONE, WORK, LDU2 ) 00267 * 00268 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) . 00269 * 00270 RESID = DLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK ) 00271 RESULT( 6 ) = ( RESID / DBLE(MAX(1,M-P)) ) / ULP 00272 * 00273 * Compute I - V1T*V1T' 00274 * 00275 CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T ) 00276 CALL DSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE, 00277 $ WORK, LDV1T ) 00278 * 00279 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) . 00280 * 00281 RESID = DLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK ) 00282 RESULT( 7 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP 00283 * 00284 * Compute I - V2T*V2T' 00285 * 00286 CALL DLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T ) 00287 CALL DSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T, 00288 $ ONE, WORK, LDV2T ) 00289 * 00290 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) . 00291 * 00292 RESID = DLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK ) 00293 RESULT( 8 ) = ( RESID / DBLE(MAX(1,M-Q)) ) / ULP 00294 * 00295 * Check sorting 00296 * 00297 RESULT(9) = REALZERO 00298 DO I = 1, R 00299 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN 00300 RESULT(9) = ULPINV 00301 END IF 00302 IF( I.GT.1) THEN 00303 IF ( THETA(I).LT.THETA(I-1) ) THEN 00304 RESULT(9) = ULPINV 00305 END IF 00306 END IF 00307 END DO 00308 * 00309 RETURN 00310 * 00311 * End of DCSDTS 00312 * 00313 END 00314
1.7.3 
