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