![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SCSDTS 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 SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, 00012 * LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, 00013 * RWORK, RESULT ) 00014 * 00015 * .. Scalar Arguments .. 00016 * INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IWORK( * ) 00020 * REAL RESULT( 9 ), RWORK( * ), THETA( * ) 00021 * REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00022 * $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ), 00023 * $ XF( LDX, * ) 00024 * .. 00025 * 00026 * 00027 *> \par Purpose: 00028 * ============= 00029 *> 00030 *> \verbatim 00031 *> 00032 *> SCSDTS tests SORCSD, which, given an M-by-M partitioned orthogonal 00033 *> matrix X, 00034 *> Q M-Q 00035 *> X = [ X11 X12 ] P , 00036 *> [ X21 X22 ] M-P 00037 *> 00038 *> computes the CSD 00039 *> 00040 *> [ U1 ]**T * [ X11 X12 ] * [ V1 ] 00041 *> [ U2 ] [ X21 X22 ] [ V2 ] 00042 *> 00043 *> [ I 0 0 | 0 0 0 ] 00044 *> [ 0 C 0 | 0 -S 0 ] 00045 *> [ 0 0 0 | 0 0 -I ] 00046 *> = [---------------------] = [ D11 D12 ] . 00047 *> [ 0 0 0 | I 0 0 ] [ D21 D22 ] 00048 *> [ 0 S 0 | 0 C 0 ] 00049 *> [ 0 0 I | 0 0 0 ] 00050 *> \endverbatim 00051 * 00052 * Arguments: 00053 * ========== 00054 * 00055 *> \param[in] M 00056 *> \verbatim 00057 *> M is INTEGER 00058 *> The number of rows of the matrix X. M >= 0. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] P 00062 *> \verbatim 00063 *> P is INTEGER 00064 *> The number of rows of the matrix X11. P >= 0. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] Q 00068 *> \verbatim 00069 *> Q is INTEGER 00070 *> The number of columns of the matrix X11. Q >= 0. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] X 00074 *> \verbatim 00075 *> X is REAL array, dimension (LDX,M) 00076 *> The M-by-M matrix X. 00077 *> \endverbatim 00078 *> 00079 *> \param[out] XF 00080 *> \verbatim 00081 *> XF is REAL array, dimension (LDX,M) 00082 *> Details of the CSD of X, as returned by SORCSD; 00083 *> see SORCSD for further details. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDX 00087 *> \verbatim 00088 *> LDX is INTEGER 00089 *> The leading dimension of the arrays X and XF. 00090 *> LDX >= max( 1,M ). 00091 *> \endverbatim 00092 *> 00093 *> \param[out] U1 00094 *> \verbatim 00095 *> U1 is REAL array, dimension(LDU1,P) 00096 *> The P-by-P orthogonal matrix U1. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] LDU1 00100 *> \verbatim 00101 *> LDU1 is INTEGER 00102 *> The leading dimension of the array U1. LDU >= max(1,P). 00103 *> \endverbatim 00104 *> 00105 *> \param[out] U2 00106 *> \verbatim 00107 *> U2 is REAL array, dimension(LDU2,M-P) 00108 *> The (M-P)-by-(M-P) orthogonal matrix U2. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDU2 00112 *> \verbatim 00113 *> LDU2 is INTEGER 00114 *> The leading dimension of the array U2. LDU >= max(1,M-P). 00115 *> \endverbatim 00116 *> 00117 *> \param[out] V1T 00118 *> \verbatim 00119 *> V1T is REAL array, dimension(LDV1T,Q) 00120 *> The Q-by-Q orthogonal matrix V1T. 00121 *> \endverbatim 00122 *> 00123 *> \param[in] LDV1T 00124 *> \verbatim 00125 *> LDV1T is INTEGER 00126 *> The leading dimension of the array V1T. LDV1T >= 00127 *> max(1,Q). 00128 *> \endverbatim 00129 *> 00130 *> \param[out] V2T 00131 *> \verbatim 00132 *> V2T is REAL array, dimension(LDV2T,M-Q) 00133 *> The (M-Q)-by-(M-Q) orthogonal matrix V2T. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] LDV2T 00137 *> \verbatim 00138 *> LDV2T is INTEGER 00139 *> The leading dimension of the array V2T. LDV2T >= 00140 *> max(1,M-Q). 00141 *> \endverbatim 00142 *> 00143 *> \param[out] THETA 00144 *> \verbatim 00145 *> THETA is REAL array, dimension MIN(P,M-P,Q,M-Q) 00146 *> The CS values of X; the essentially diagonal matrices C and 00147 *> S are constructed from THETA; see subroutine SORCSD for 00148 *> details. 00149 *> \endverbatim 00150 *> 00151 *> \param[out] IWORK 00152 *> \verbatim 00153 *> IWORK is INTEGER array, dimension (M) 00154 *> \endverbatim 00155 *> 00156 *> \param[out] WORK 00157 *> \verbatim 00158 *> WORK is REAL array, dimension (LWORK) 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LWORK 00162 *> \verbatim 00163 *> LWORK is INTEGER 00164 *> The dimension of the array WORK 00165 *> \endverbatim 00166 *> 00167 *> \param[out] RWORK 00168 *> \verbatim 00169 *> RWORK is REAL array 00170 *> \endverbatim 00171 *> 00172 *> \param[out] RESULT 00173 *> \verbatim 00174 *> RESULT is REAL array, dimension (9) 00175 *> The test ratios: 00176 *> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) 00177 *> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) 00178 *> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) 00179 *> RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) 00180 *> RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP ) 00181 *> RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP ) 00182 *> RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP ) 00183 *> RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP ) 00184 *> RESULT(9) = 0 if THETA is in increasing order and 00185 *> all angles are in [0,pi/2]; 00186 *> = ULPINV otherwise. 00187 *> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). ) 00188 *> \endverbatim 00189 * 00190 * Authors: 00191 * ======== 00192 * 00193 *> \author Univ. of Tennessee 00194 *> \author Univ. of California Berkeley 00195 *> \author Univ. of Colorado Denver 00196 *> \author NAG Ltd. 00197 * 00198 *> \date November 2011 00199 * 00200 *> \ingroup single_eig 00201 * 00202 * ===================================================================== 00203 SUBROUTINE SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, 00204 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, 00205 $ RWORK, RESULT ) 00206 * 00207 * -- LAPACK test routine (version 3.4.0) -- 00208 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00210 * November 2011 00211 * 00212 * .. Scalar Arguments .. 00213 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q 00214 * .. 00215 * .. Array Arguments .. 00216 INTEGER IWORK( * ) 00217 REAL RESULT( 9 ), RWORK( * ), THETA( * ) 00218 REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 00219 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ), 00220 $ XF( LDX, * ) 00221 * .. 00222 * 00223 * ===================================================================== 00224 * 00225 * .. Parameters .. 00226 REAL PIOVER2, REALONE, REALZERO 00227 PARAMETER ( PIOVER2 = 1.57079632679489662E0, 00228 $ REALONE = 1.0E0, REALZERO = 0.0E0 ) 00229 REAL ZERO, ONE 00230 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00231 * .. 00232 * .. Local Scalars .. 00233 INTEGER I, INFO, R 00234 REAL EPS2, RESID, ULP, ULPINV 00235 * .. 00236 * .. External Functions .. 00237 REAL SLAMCH, SLANGE, SLANSY 00238 EXTERNAL SLAMCH, SLANGE, SLANSY 00239 * .. 00240 * .. External Subroutines .. 00241 EXTERNAL SGEMM, SLACPY, SLASET, SORCSD, SSYRK 00242 * .. 00243 * .. Intrinsic Functions .. 00244 INTRINSIC REAL, MAX, MIN 00245 * .. 00246 * .. Executable Statements .. 00247 * 00248 ULP = SLAMCH( 'Precision' ) 00249 ULPINV = REALONE / ULP 00250 CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX ) 00251 CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX, 00252 $ ONE, WORK, LDX ) 00253 IF (M.GT.0) THEN 00254 EPS2 = MAX( ULP, 00255 $ SLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) ) 00256 ELSE 00257 EPS2 = ULP 00258 END IF 00259 R = MIN( P, M-P, Q, M-Q ) 00260 * 00261 * Copy the matrix X to the array XF. 00262 * 00263 CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX ) 00264 * 00265 * Compute the CSD 00266 * 00267 CALL SORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX, 00268 $ XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX, 00269 $ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 00270 $ WORK, LWORK, IWORK, INFO ) 00271 * 00272 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22] 00273 * 00274 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE, 00275 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00276 * 00277 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE, 00278 $ U1, LDU1, WORK, LDX, ZERO, X, LDX ) 00279 * 00280 DO I = 1, MIN(P,Q)-R 00281 X(I,I) = X(I,I) - ONE 00282 END DO 00283 DO I = 1, R 00284 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = 00285 $ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I)) 00286 END DO 00287 * 00288 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q, 00289 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00290 * 00291 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P, 00292 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX ) 00293 * 00294 DO I = 1, MIN(P,M-Q)-R 00295 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE 00296 END DO 00297 DO I = 1, R 00298 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) = 00299 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) + 00300 $ SIN(THETA(R-I+1)) 00301 END DO 00302 * 00303 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE, 00304 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX ) 00305 * 00306 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P, 00307 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX ) 00308 * 00309 DO I = 1, MIN(M-P,Q)-R 00310 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE 00311 END DO 00312 DO I = 1, R 00313 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) = 00314 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) - 00315 $ SIN(THETA(R-I+1)) 00316 END DO 00317 * 00318 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q, 00319 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX ) 00320 * 00321 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P, 00322 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX ) 00323 * 00324 DO I = 1, MIN(M-P,M-Q)-R 00325 X(P+I,Q+I) = X(P+I,Q+I) - ONE 00326 END DO 00327 DO I = 1, R 00328 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) = 00329 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) - 00330 $ COS(THETA(I)) 00331 END DO 00332 * 00333 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) . 00334 * 00335 RESID = SLANGE( '1', P, Q, X, LDX, RWORK ) 00336 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2 00337 * 00338 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) . 00339 * 00340 RESID = SLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK ) 00341 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2 00342 * 00343 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) . 00344 * 00345 RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK ) 00346 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2 00347 * 00348 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) . 00349 * 00350 RESID = SLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK ) 00351 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2 00352 * 00353 * Compute I - U1'*U1 00354 * 00355 CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 ) 00356 CALL SSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1, 00357 $ ONE, WORK, LDU1 ) 00358 * 00359 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) . 00360 * 00361 RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK ) 00362 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP 00363 * 00364 * Compute I - U2'*U2 00365 * 00366 CALL SLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 ) 00367 CALL SSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2, 00368 $ LDU2, ONE, WORK, LDU2 ) 00369 * 00370 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) . 00371 * 00372 RESID = SLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK ) 00373 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP 00374 * 00375 * Compute I - V1T*V1T' 00376 * 00377 CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T ) 00378 CALL SSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE, 00379 $ WORK, LDV1T ) 00380 * 00381 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) . 00382 * 00383 RESID = SLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK ) 00384 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP 00385 * 00386 * Compute I - V2T*V2T' 00387 * 00388 CALL SLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T ) 00389 CALL SSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T, 00390 $ ONE, WORK, LDV2T ) 00391 * 00392 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) . 00393 * 00394 RESID = SLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK ) 00395 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP 00396 * 00397 * Check sorting 00398 * 00399 RESULT(9) = REALZERO 00400 DO I = 1, R 00401 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN 00402 RESULT(9) = ULPINV 00403 END IF 00404 IF( I.GT.1) THEN 00405 IF ( THETA(I).LT.THETA(I-1) ) THEN 00406 RESULT(9) = ULPINV 00407 END IF 00408 END IF 00409 END DO 00410 * 00411 RETURN 00412 * 00413 * End of SCSDTS 00414 * 00415 END 00416