|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE SGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, 00002 $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, 00003 $ LWORK, RWORK, RESULT ) 00004 * 00005 * -- LAPACK test routine (version 3.1) -- 00006 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IWORK( * ) 00014 REAL A( LDA, * ), AF( LDA, * ), ALPHA( * ), 00015 $ B( LDB, * ), BETA( * ), BF( LDB, * ), 00016 $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), 00017 $ RWORK( * ), U( LDU, * ), V( LDV, * ), 00018 $ WORK( LWORK ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * SGSVTS tests SGGSVD, which computes the GSVD of an M-by-N matrix A 00025 * and a P-by-N matrix B: 00026 * U'*A*Q = D1*R and V'*B*Q = D2*R. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * M (input) INTEGER 00032 * The number of rows of the matrix A. M >= 0. 00033 * 00034 * P (input) INTEGER 00035 * The number of rows of the matrix B. P >= 0. 00036 * 00037 * N (input) INTEGER 00038 * The number of columns of the matrices A and B. N >= 0. 00039 * 00040 * A (input) REAL array, dimension (LDA,M) 00041 * The M-by-N matrix A. 00042 * 00043 * AF (output) REAL array, dimension (LDA,N) 00044 * Details of the GSVD of A and B, as returned by SGGSVD, 00045 * see SGGSVD for further details. 00046 * 00047 * LDA (input) INTEGER 00048 * The leading dimension of the arrays A and AF. 00049 * LDA >= max( 1,M ). 00050 * 00051 * B (input) REAL array, dimension (LDB,P) 00052 * On entry, the P-by-N matrix B. 00053 * 00054 * BF (output) REAL array, dimension (LDB,N) 00055 * Details of the GSVD of A and B, as returned by SGGSVD, 00056 * see SGGSVD for further details. 00057 * 00058 * LDB (input) INTEGER 00059 * The leading dimension of the arrays B and BF. 00060 * LDB >= max(1,P). 00061 * 00062 * U (output) REAL array, dimension(LDU,M) 00063 * The M by M orthogonal matrix U. 00064 * 00065 * LDU (input) INTEGER 00066 * The leading dimension of the array U. LDU >= max(1,M). 00067 * 00068 * V (output) REAL array, dimension(LDV,M) 00069 * The P by P orthogonal matrix V. 00070 * 00071 * LDV (input) INTEGER 00072 * The leading dimension of the array V. LDV >= max(1,P). 00073 * 00074 * Q (output) REAL array, dimension(LDQ,N) 00075 * The N by N orthogonal matrix Q. 00076 * 00077 * LDQ (input) INTEGER 00078 * The leading dimension of the array Q. LDQ >= max(1,N). 00079 * 00080 * ALPHA (output) REAL array, dimension (N) 00081 * BETA (output) REAL array, dimension (N) 00082 * The generalized singular value pairs of A and B, the 00083 * ``diagonal'' matrices D1 and D2 are constructed from 00084 * ALPHA and BETA, see subroutine SGGSVD for details. 00085 * 00086 * R (output) REAL array, dimension(LDQ,N) 00087 * The upper triangular matrix R. 00088 * 00089 * LDR (input) INTEGER 00090 * The leading dimension of the array R. LDR >= max(1,N). 00091 * 00092 * IWORK (workspace) INTEGER array, dimension (N) 00093 * 00094 * WORK (workspace) REAL array, dimension (LWORK) 00095 * 00096 * LWORK (input) INTEGER 00097 * The dimension of the array WORK, 00098 * LWORK >= max(M,P,N)*max(M,P,N). 00099 * 00100 * RWORK (workspace) REAL array, dimension (max(M,P,N)) 00101 * 00102 * RESULT (output) REAL array, dimension (6) 00103 * The test ratios: 00104 * RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) 00105 * RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) 00106 * RESULT(3) = norm( I - U'*U ) / ( M*ULP ) 00107 * RESULT(4) = norm( I - V'*V ) / ( P*ULP ) 00108 * RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) 00109 * RESULT(6) = 0 if ALPHA is in decreasing order; 00110 * = ULPINV otherwise. 00111 * 00112 * ===================================================================== 00113 * 00114 * .. Parameters .. 00115 REAL ZERO, ONE 00116 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00117 * .. 00118 * .. Local Scalars .. 00119 INTEGER I, INFO, J, K, L 00120 REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL 00121 * .. 00122 * .. External Functions .. 00123 REAL SLAMCH, SLANGE, SLANSY 00124 EXTERNAL SLAMCH, SLANGE, SLANSY 00125 * .. 00126 * .. External Subroutines .. 00127 EXTERNAL SCOPY, SGEMM, SGGSVD, SLACPY, SLASET, SSYRK 00128 * .. 00129 * .. Intrinsic Functions .. 00130 INTRINSIC MAX, MIN, REAL 00131 * .. 00132 * .. Executable Statements .. 00133 * 00134 ULP = SLAMCH( 'Precision' ) 00135 ULPINV = ONE / ULP 00136 UNFL = SLAMCH( 'Safe minimum' ) 00137 * 00138 * Copy the matrix A to the array AF. 00139 * 00140 CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA ) 00141 CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB ) 00142 * 00143 ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 00144 BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 00145 * 00146 * Factorize the matrices A and B in the arrays AF and BF. 00147 * 00148 CALL SGGSVD( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, 00149 $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, 00150 $ INFO ) 00151 * 00152 * Copy R 00153 * 00154 DO 20 I = 1, MIN( K+L, M ) 00155 DO 10 J = I, K + L 00156 R( I, J ) = AF( I, N-K-L+J ) 00157 10 CONTINUE 00158 20 CONTINUE 00159 * 00160 IF( M-K-L.LT.0 ) THEN 00161 DO 40 I = M + 1, K + L 00162 DO 30 J = I, K + L 00163 R( I, J ) = BF( I-K, N-K-L+J ) 00164 30 CONTINUE 00165 40 CONTINUE 00166 END IF 00167 * 00168 * Compute A:= U'*A*Q - D1*R 00169 * 00170 CALL SGEMM( 'No transpose', 'No transpose', M, N, N, ONE, A, LDA, 00171 $ Q, LDQ, ZERO, WORK, LDA ) 00172 * 00173 CALL SGEMM( 'Transpose', 'No transpose', M, N, M, ONE, U, LDU, 00174 $ WORK, LDA, ZERO, A, LDA ) 00175 * 00176 DO 60 I = 1, K 00177 DO 50 J = I, K + L 00178 A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) 00179 50 CONTINUE 00180 60 CONTINUE 00181 * 00182 DO 80 I = K + 1, MIN( K+L, M ) 00183 DO 70 J = I, K + L 00184 A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) 00185 70 CONTINUE 00186 80 CONTINUE 00187 * 00188 * Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . 00189 * 00190 RESID = SLANGE( '1', M, N, A, LDA, RWORK ) 00191 * 00192 IF( ANORM.GT.ZERO ) THEN 00193 RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) / 00194 $ ULP 00195 ELSE 00196 RESULT( 1 ) = ZERO 00197 END IF 00198 * 00199 * Compute B := V'*B*Q - D2*R 00200 * 00201 CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, B, LDB, 00202 $ Q, LDQ, ZERO, WORK, LDB ) 00203 * 00204 CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, V, LDV, 00205 $ WORK, LDB, ZERO, B, LDB ) 00206 * 00207 DO 100 I = 1, L 00208 DO 90 J = I, L 00209 B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) 00210 90 CONTINUE 00211 100 CONTINUE 00212 * 00213 * Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . 00214 * 00215 RESID = SLANGE( '1', P, N, B, LDB, RWORK ) 00216 IF( BNORM.GT.ZERO ) THEN 00217 RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) / 00218 $ ULP 00219 ELSE 00220 RESULT( 2 ) = ZERO 00221 END IF 00222 * 00223 * Compute I - U'*U 00224 * 00225 CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDQ ) 00226 CALL SSYRK( 'Upper', 'Transpose', M, M, -ONE, U, LDU, ONE, WORK, 00227 $ LDU ) 00228 * 00229 * Compute norm( I - U'*U ) / ( M * ULP ) . 00230 * 00231 RESID = SLANSY( '1', 'Upper', M, WORK, LDU, RWORK ) 00232 RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP 00233 * 00234 * Compute I - V'*V 00235 * 00236 CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDV ) 00237 CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, V, LDV, ONE, WORK, 00238 $ LDV ) 00239 * 00240 * Compute norm( I - V'*V ) / ( P * ULP ) . 00241 * 00242 RESID = SLANSY( '1', 'Upper', P, WORK, LDV, RWORK ) 00243 RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP 00244 * 00245 * Compute I - Q'*Q 00246 * 00247 CALL SLASET( 'Full', N, N, ZERO, ONE, WORK, LDQ ) 00248 CALL SSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDQ, ONE, WORK, 00249 $ LDQ ) 00250 * 00251 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00252 * 00253 RESID = SLANSY( '1', 'Upper', N, WORK, LDQ, RWORK ) 00254 RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP 00255 * 00256 * Check sorting 00257 * 00258 CALL SCOPY( N, ALPHA, 1, WORK, 1 ) 00259 DO 110 I = K + 1, MIN( K+L, M ) 00260 J = IWORK( I ) 00261 IF( I.NE.J ) THEN 00262 TEMP = WORK( I ) 00263 WORK( I ) = WORK( J ) 00264 WORK( J ) = TEMP 00265 END IF 00266 110 CONTINUE 00267 * 00268 RESULT( 6 ) = ZERO 00269 DO 120 I = K + 1, MIN( K+L, M ) - 1 00270 IF( WORK( I ).LT.WORK( I+1 ) ) 00271 $ RESULT( 6 ) = ULPINV 00272 120 CONTINUE 00273 * 00274 RETURN 00275 * 00276 * End of SGSVTS 00277 * 00278 END
1.7.3 
