1*> \brief \b SGSVTS3 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 12* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 13* LWORK, RWORK, RESULT ) 14* 15* .. Scalar Arguments .. 16* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 17* .. 18* .. Array Arguments .. 19* INTEGER IWORK( * ) 20* REAL A( LDA, * ), AF( LDA, * ), ALPHA( * ), 21* $ B( LDB, * ), BETA( * ), BF( LDB, * ), 22* $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), 23* $ RWORK( * ), U( LDU, * ), V( LDV, * ), 24* $ WORK( LWORK ) 25* .. 26* 27* 28*> \par Purpose: 29* ============= 30*> 31*> \verbatim 32*> 33*> SGSVTS3 tests SGGSVD3, which computes the GSVD of an M-by-N matrix A 34*> and a P-by-N matrix B: 35*> U'*A*Q = D1*R and V'*B*Q = D2*R. 36*> \endverbatim 37* 38* Arguments: 39* ========== 40* 41*> \param[in] M 42*> \verbatim 43*> M is INTEGER 44*> The number of rows of the matrix A. M >= 0. 45*> \endverbatim 46*> 47*> \param[in] P 48*> \verbatim 49*> P is INTEGER 50*> The number of rows of the matrix B. P >= 0. 51*> \endverbatim 52*> 53*> \param[in] N 54*> \verbatim 55*> N is INTEGER 56*> The number of columns of the matrices A and B. N >= 0. 57*> \endverbatim 58*> 59*> \param[in] A 60*> \verbatim 61*> A is REAL array, dimension (LDA,M) 62*> The M-by-N matrix A. 63*> \endverbatim 64*> 65*> \param[out] AF 66*> \verbatim 67*> AF is REAL array, dimension (LDA,N) 68*> Details of the GSVD of A and B, as returned by SGGSVD3, 69*> see SGGSVD3 for further details. 70*> \endverbatim 71*> 72*> \param[in] LDA 73*> \verbatim 74*> LDA is INTEGER 75*> The leading dimension of the arrays A and AF. 76*> LDA >= max( 1,M ). 77*> \endverbatim 78*> 79*> \param[in] B 80*> \verbatim 81*> B is REAL array, dimension (LDB,P) 82*> On entry, the P-by-N matrix B. 83*> \endverbatim 84*> 85*> \param[out] BF 86*> \verbatim 87*> BF is REAL array, dimension (LDB,N) 88*> Details of the GSVD of A and B, as returned by SGGSVD3, 89*> see SGGSVD3 for further details. 90*> \endverbatim 91*> 92*> \param[in] LDB 93*> \verbatim 94*> LDB is INTEGER 95*> The leading dimension of the arrays B and BF. 96*> LDB >= max(1,P). 97*> \endverbatim 98*> 99*> \param[out] U 100*> \verbatim 101*> U is REAL array, dimension(LDU,M) 102*> The M by M orthogonal matrix U. 103*> \endverbatim 104*> 105*> \param[in] LDU 106*> \verbatim 107*> LDU is INTEGER 108*> The leading dimension of the array U. LDU >= max(1,M). 109*> \endverbatim 110*> 111*> \param[out] V 112*> \verbatim 113*> V is REAL array, dimension(LDV,M) 114*> The P by P orthogonal matrix V. 115*> \endverbatim 116*> 117*> \param[in] LDV 118*> \verbatim 119*> LDV is INTEGER 120*> The leading dimension of the array V. LDV >= max(1,P). 121*> \endverbatim 122*> 123*> \param[out] Q 124*> \verbatim 125*> Q is REAL array, dimension(LDQ,N) 126*> The N by N orthogonal matrix Q. 127*> \endverbatim 128*> 129*> \param[in] LDQ 130*> \verbatim 131*> LDQ is INTEGER 132*> The leading dimension of the array Q. LDQ >= max(1,N). 133*> \endverbatim 134*> 135*> \param[out] ALPHA 136*> \verbatim 137*> ALPHA is REAL array, dimension (N) 138*> \endverbatim 139*> 140*> \param[out] BETA 141*> \verbatim 142*> BETA is REAL array, dimension (N) 143*> 144*> The generalized singular value pairs of A and B, the 145*> ``diagonal'' matrices D1 and D2 are constructed from 146*> ALPHA and BETA, see subroutine SGGSVD3 for details. 147*> \endverbatim 148*> 149*> \param[out] R 150*> \verbatim 151*> R is REAL array, dimension(LDQ,N) 152*> The upper triangular matrix R. 153*> \endverbatim 154*> 155*> \param[in] LDR 156*> \verbatim 157*> LDR is INTEGER 158*> The leading dimension of the array R. LDR >= max(1,N). 159*> \endverbatim 160*> 161*> \param[out] IWORK 162*> \verbatim 163*> IWORK is INTEGER array, dimension (N) 164*> \endverbatim 165*> 166*> \param[out] WORK 167*> \verbatim 168*> WORK is REAL array, dimension (LWORK) 169*> \endverbatim 170*> 171*> \param[in] LWORK 172*> \verbatim 173*> LWORK is INTEGER 174*> The dimension of the array WORK, 175*> LWORK >= max(M,P,N)*max(M,P,N). 176*> \endverbatim 177*> 178*> \param[out] RWORK 179*> \verbatim 180*> RWORK is REAL array, dimension (max(M,P,N)) 181*> \endverbatim 182*> 183*> \param[out] RESULT 184*> \verbatim 185*> RESULT is REAL array, dimension (6) 186*> The test ratios: 187*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) 188*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) 189*> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) 190*> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) 191*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) 192*> RESULT(6) = 0 if ALPHA is in decreasing order; 193*> = ULPINV otherwise. 194*> \endverbatim 195* 196* Authors: 197* ======== 198* 199*> \author Univ. of Tennessee 200*> \author Univ. of California Berkeley 201*> \author Univ. of Colorado Denver 202*> \author NAG Ltd. 203* 204*> \ingroup single_eig 205* 206* ===================================================================== 207 SUBROUTINE SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 208 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 209 $ LWORK, RWORK, RESULT ) 210* 211* -- LAPACK test routine -- 212* -- LAPACK is a software package provided by Univ. of Tennessee, -- 213* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 214* 215* .. Scalar Arguments .. 216 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 217* .. 218* .. Array Arguments .. 219 INTEGER IWORK( * ) 220 REAL A( LDA, * ), AF( LDA, * ), ALPHA( * ), 221 $ B( LDB, * ), BETA( * ), BF( LDB, * ), 222 $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), 223 $ RWORK( * ), U( LDU, * ), V( LDV, * ), 224 $ WORK( LWORK ) 225* .. 226* 227* ===================================================================== 228* 229* .. Parameters .. 230 REAL ZERO, ONE 231 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 232* .. 233* .. Local Scalars .. 234 INTEGER I, INFO, J, K, L 235 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL 236* .. 237* .. External Functions .. 238 REAL SLAMCH, SLANGE, SLANSY 239 EXTERNAL SLAMCH, SLANGE, SLANSY 240* .. 241* .. External Subroutines .. 242 EXTERNAL SCOPY, SGEMM, SGGSVD3, SLACPY, SLASET, SSYRK 243* .. 244* .. Intrinsic Functions .. 245 INTRINSIC MAX, MIN, REAL 246* .. 247* .. Executable Statements .. 248* 249 ULP = SLAMCH( 'Precision' ) 250 ULPINV = ONE / ULP 251 UNFL = SLAMCH( 'Safe minimum' ) 252* 253* Copy the matrix A to the array AF. 254* 255 CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA ) 256 CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB ) 257* 258 ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 259 BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 260* 261* Factorize the matrices A and B in the arrays AF and BF. 262* 263 CALL SGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, 264 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, 265 $ IWORK, INFO ) 266* 267* Copy R 268* 269 DO 20 I = 1, MIN( K+L, M ) 270 DO 10 J = I, K + L 271 R( I, J ) = AF( I, N-K-L+J ) 272 10 CONTINUE 273 20 CONTINUE 274* 275 IF( M-K-L.LT.0 ) THEN 276 DO 40 I = M + 1, K + L 277 DO 30 J = I, K + L 278 R( I, J ) = BF( I-K, N-K-L+J ) 279 30 CONTINUE 280 40 CONTINUE 281 END IF 282* 283* Compute A:= U'*A*Q - D1*R 284* 285 CALL SGEMM( 'No transpose', 'No transpose', M, N, N, ONE, A, LDA, 286 $ Q, LDQ, ZERO, WORK, LDA ) 287* 288 CALL SGEMM( 'Transpose', 'No transpose', M, N, M, ONE, U, LDU, 289 $ WORK, LDA, ZERO, A, LDA ) 290* 291 DO 60 I = 1, K 292 DO 50 J = I, K + L 293 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) 294 50 CONTINUE 295 60 CONTINUE 296* 297 DO 80 I = K + 1, MIN( K+L, M ) 298 DO 70 J = I, K + L 299 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) 300 70 CONTINUE 301 80 CONTINUE 302* 303* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . 304* 305 RESID = SLANGE( '1', M, N, A, LDA, RWORK ) 306* 307 IF( ANORM.GT.ZERO ) THEN 308 RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) / 309 $ ULP 310 ELSE 311 RESULT( 1 ) = ZERO 312 END IF 313* 314* Compute B := V'*B*Q - D2*R 315* 316 CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, B, LDB, 317 $ Q, LDQ, ZERO, WORK, LDB ) 318* 319 CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, V, LDV, 320 $ WORK, LDB, ZERO, B, LDB ) 321* 322 DO 100 I = 1, L 323 DO 90 J = I, L 324 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) 325 90 CONTINUE 326 100 CONTINUE 327* 328* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . 329* 330 RESID = SLANGE( '1', P, N, B, LDB, RWORK ) 331 IF( BNORM.GT.ZERO ) THEN 332 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) / 333 $ ULP 334 ELSE 335 RESULT( 2 ) = ZERO 336 END IF 337* 338* Compute I - U'*U 339* 340 CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDQ ) 341 CALL SSYRK( 'Upper', 'Transpose', M, M, -ONE, U, LDU, ONE, WORK, 342 $ LDU ) 343* 344* Compute norm( I - U'*U ) / ( M * ULP ) . 345* 346 RESID = SLANSY( '1', 'Upper', M, WORK, LDU, RWORK ) 347 RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP 348* 349* Compute I - V'*V 350* 351 CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDV ) 352 CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, V, LDV, ONE, WORK, 353 $ LDV ) 354* 355* Compute norm( I - V'*V ) / ( P * ULP ) . 356* 357 RESID = SLANSY( '1', 'Upper', P, WORK, LDV, RWORK ) 358 RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP 359* 360* Compute I - Q'*Q 361* 362 CALL SLASET( 'Full', N, N, ZERO, ONE, WORK, LDQ ) 363 CALL SSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDQ, ONE, WORK, 364 $ LDQ ) 365* 366* Compute norm( I - Q'*Q ) / ( N * ULP ) . 367* 368 RESID = SLANSY( '1', 'Upper', N, WORK, LDQ, RWORK ) 369 RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP 370* 371* Check sorting 372* 373 CALL SCOPY( N, ALPHA, 1, WORK, 1 ) 374 DO 110 I = K + 1, MIN( K+L, M ) 375 J = IWORK( I ) 376 IF( I.NE.J ) THEN 377 TEMP = WORK( I ) 378 WORK( I ) = WORK( J ) 379 WORK( J ) = TEMP 380 END IF 381 110 CONTINUE 382* 383 RESULT( 6 ) = ZERO 384 DO 120 I = K + 1, MIN( K+L, M ) - 1 385 IF( WORK( I ).LT.WORK( I+1 ) ) 386 $ RESULT( 6 ) = ULPINV 387 120 CONTINUE 388* 389 RETURN 390* 391* End of SGSVTS3 392* 393 END 394