![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGQRTS 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 SGQRTS( N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, 00012 * BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER LDA, LDB, LWORK, M, P, N 00016 * .. 00017 * .. Array Arguments .. 00018 * REAL A( LDA, * ), AF( LDA, * ), R( LDA, * ), 00019 * $ Q( LDA, * ), B( LDB, * ), BF( LDB, * ), 00020 * $ T( LDB, * ), Z( LDB, * ), BWK( LDB, * ), 00021 * $ TAUA( * ), TAUB( * ), RESULT( 4 ), 00022 * $ RWORK( * ), WORK( LWORK ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> SGQRTS tests SGGQRF, which computes the GQR factorization of an 00032 *> N-by-M matrix A and a N-by-P matrix B: A = Q*R and B = Q*T*Z. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] N 00039 *> \verbatim 00040 *> N is INTEGER 00041 *> The number of rows of the matrices A and B. N >= 0. 00042 *> \endverbatim 00043 *> 00044 *> \param[in] M 00045 *> \verbatim 00046 *> M is INTEGER 00047 *> The number of columns of the matrix A. M >= 0. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] P 00051 *> \verbatim 00052 *> P is INTEGER 00053 *> The number of columns of the matrix B. P >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] A 00057 *> \verbatim 00058 *> A is REAL array, dimension (LDA,M) 00059 *> The N-by-M matrix A. 00060 *> \endverbatim 00061 *> 00062 *> \param[out] AF 00063 *> \verbatim 00064 *> AF is REAL array, dimension (LDA,N) 00065 *> Details of the GQR factorization of A and B, as returned 00066 *> by SGGQRF, see SGGQRF for further details. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] Q 00070 *> \verbatim 00071 *> Q is REAL array, dimension (LDA,N) 00072 *> The M-by-M orthogonal matrix Q. 00073 *> \endverbatim 00074 *> 00075 *> \param[out] R 00076 *> \verbatim 00077 *> R is REAL array, dimension (LDA,MAX(M,N)) 00078 *> \endverbatim 00079 *> 00080 *> \param[in] LDA 00081 *> \verbatim 00082 *> LDA is INTEGER 00083 *> The leading dimension of the arrays A, AF, R and Q. 00084 *> LDA >= max(M,N). 00085 *> \endverbatim 00086 *> 00087 *> \param[out] TAUA 00088 *> \verbatim 00089 *> TAUA is REAL array, dimension (min(M,N)) 00090 *> The scalar factors of the elementary reflectors, as returned 00091 *> by SGGQRF. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] B 00095 *> \verbatim 00096 *> B is REAL array, dimension (LDB,P) 00097 *> On entry, the N-by-P matrix A. 00098 *> \endverbatim 00099 *> 00100 *> \param[out] BF 00101 *> \verbatim 00102 *> BF is REAL array, dimension (LDB,N) 00103 *> Details of the GQR factorization of A and B, as returned 00104 *> by SGGQRF, see SGGQRF for further details. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Z 00108 *> \verbatim 00109 *> Z is REAL array, dimension (LDB,P) 00110 *> The P-by-P orthogonal matrix Z. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] T 00114 *> \verbatim 00115 *> T is REAL array, dimension (LDB,max(P,N)) 00116 *> \endverbatim 00117 *> 00118 *> \param[out] BWK 00119 *> \verbatim 00120 *> BWK is REAL array, dimension (LDB,N) 00121 *> \endverbatim 00122 *> 00123 *> \param[in] LDB 00124 *> \verbatim 00125 *> LDB is INTEGER 00126 *> The leading dimension of the arrays B, BF, Z and T. 00127 *> LDB >= max(P,N). 00128 *> \endverbatim 00129 *> 00130 *> \param[out] TAUB 00131 *> \verbatim 00132 *> TAUB is REAL array, dimension (min(P,N)) 00133 *> The scalar factors of the elementary reflectors, as returned 00134 *> by SGGRQF. 00135 *> \endverbatim 00136 *> 00137 *> \param[out] WORK 00138 *> \verbatim 00139 *> WORK is REAL array, dimension (LWORK) 00140 *> \endverbatim 00141 *> 00142 *> \param[in] LWORK 00143 *> \verbatim 00144 *> LWORK is INTEGER 00145 *> The dimension of the array WORK, LWORK >= max(N,M,P)**2. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] RWORK 00149 *> \verbatim 00150 *> RWORK is REAL array, dimension (max(N,M,P)) 00151 *> \endverbatim 00152 *> 00153 *> \param[out] RESULT 00154 *> \verbatim 00155 *> RESULT is REAL array, dimension (4) 00156 *> The test ratios: 00157 *> RESULT(1) = norm( R - Q'*A ) / ( MAX(M,N)*norm(A)*ULP) 00158 *> RESULT(2) = norm( T*Z - Q'*B ) / (MAX(P,N)*norm(B)*ULP) 00159 *> RESULT(3) = norm( I - Q'*Q ) / ( M*ULP ) 00160 *> RESULT(4) = norm( I - Z'*Z ) / ( P*ULP ) 00161 *> \endverbatim 00162 * 00163 * Authors: 00164 * ======== 00165 * 00166 *> \author Univ. of Tennessee 00167 *> \author Univ. of California Berkeley 00168 *> \author Univ. of Colorado Denver 00169 *> \author NAG Ltd. 00170 * 00171 *> \date November 2011 00172 * 00173 *> \ingroup single_eig 00174 * 00175 * ===================================================================== 00176 SUBROUTINE SGQRTS( N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, 00177 $ BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT ) 00178 * 00179 * -- LAPACK test routine (version 3.4.0) -- 00180 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00181 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00182 * November 2011 00183 * 00184 * .. Scalar Arguments .. 00185 INTEGER LDA, LDB, LWORK, M, P, N 00186 * .. 00187 * .. Array Arguments .. 00188 REAL A( LDA, * ), AF( LDA, * ), R( LDA, * ), 00189 $ Q( LDA, * ), B( LDB, * ), BF( LDB, * ), 00190 $ T( LDB, * ), Z( LDB, * ), BWK( LDB, * ), 00191 $ TAUA( * ), TAUB( * ), RESULT( 4 ), 00192 $ RWORK( * ), WORK( LWORK ) 00193 * .. 00194 * 00195 * ===================================================================== 00196 * 00197 * .. Parameters .. 00198 REAL ZERO, ONE 00199 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00200 REAL ROGUE 00201 PARAMETER ( ROGUE = -1.0E+10 ) 00202 * .. 00203 * .. Local Scalars .. 00204 INTEGER INFO 00205 REAL ANORM, BNORM, ULP, UNFL, RESID 00206 * .. 00207 * .. External Functions .. 00208 REAL SLAMCH, SLANGE, SLANSY 00209 EXTERNAL SLAMCH, SLANGE, SLANSY 00210 * .. 00211 * .. External Subroutines .. 00212 EXTERNAL SGEMM, SLACPY, SLASET, SORGQR, 00213 $ SORGRQ, SSYRK 00214 * .. 00215 * .. Intrinsic Functions .. 00216 INTRINSIC MAX, MIN, REAL 00217 * .. 00218 * .. Executable Statements .. 00219 * 00220 ULP = SLAMCH( 'Precision' ) 00221 UNFL = SLAMCH( 'Safe minimum' ) 00222 * 00223 * Copy the matrix A to the array AF. 00224 * 00225 CALL SLACPY( 'Full', N, M, A, LDA, AF, LDA ) 00226 CALL SLACPY( 'Full', N, P, B, LDB, BF, LDB ) 00227 * 00228 ANORM = MAX( SLANGE( '1', N, M, A, LDA, RWORK ), UNFL ) 00229 BNORM = MAX( SLANGE( '1', N, P, B, LDB, RWORK ), UNFL ) 00230 * 00231 * Factorize the matrices A and B in the arrays AF and BF. 00232 * 00233 CALL SGGQRF( N, M, P, AF, LDA, TAUA, BF, LDB, TAUB, WORK, 00234 $ LWORK, INFO ) 00235 * 00236 * Generate the N-by-N matrix Q 00237 * 00238 CALL SLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA ) 00239 CALL SLACPY( 'Lower', N-1, M, AF( 2,1 ), LDA, Q( 2,1 ), LDA ) 00240 CALL SORGQR( N, N, MIN( N, M ), Q, LDA, TAUA, WORK, LWORK, INFO ) 00241 * 00242 * Generate the P-by-P matrix Z 00243 * 00244 CALL SLASET( 'Full', P, P, ROGUE, ROGUE, Z, LDB ) 00245 IF( N.LE.P ) THEN 00246 IF( N.GT.0 .AND. N.LT.P ) 00247 $ CALL SLACPY( 'Full', N, P-N, BF, LDB, Z( P-N+1, 1 ), LDB ) 00248 IF( N.GT.1 ) 00249 $ CALL SLACPY( 'Lower', N-1, N-1, BF( 2, P-N+1 ), LDB, 00250 $ Z( P-N+2, P-N+1 ), LDB ) 00251 ELSE 00252 IF( P.GT.1) 00253 $ CALL SLACPY( 'Lower', P-1, P-1, BF( N-P+2, 1 ), LDB, 00254 $ Z( 2, 1 ), LDB ) 00255 END IF 00256 CALL SORGRQ( P, P, MIN( N, P ), Z, LDB, TAUB, WORK, LWORK, INFO ) 00257 * 00258 * Copy R 00259 * 00260 CALL SLASET( 'Full', N, M, ZERO, ZERO, R, LDA ) 00261 CALL SLACPY( 'Upper', N, M, AF, LDA, R, LDA ) 00262 * 00263 * Copy T 00264 * 00265 CALL SLASET( 'Full', N, P, ZERO, ZERO, T, LDB ) 00266 IF( N.LE.P ) THEN 00267 CALL SLACPY( 'Upper', N, N, BF( 1, P-N+1 ), LDB, T( 1, P-N+1 ), 00268 $ LDB ) 00269 ELSE 00270 CALL SLACPY( 'Full', N-P, P, BF, LDB, T, LDB ) 00271 CALL SLACPY( 'Upper', P, P, BF( N-P+1, 1 ), LDB, T( N-P+1, 1 ), 00272 $ LDB ) 00273 END IF 00274 * 00275 * Compute R - Q'*A 00276 * 00277 CALL SGEMM( 'Transpose', 'No transpose', N, M, N, -ONE, Q, LDA, A, 00278 $ LDA, ONE, R, LDA ) 00279 * 00280 * Compute norm( R - Q'*A ) / ( MAX(M,N)*norm(A)*ULP ) . 00281 * 00282 RESID = SLANGE( '1', N, M, R, LDA, RWORK ) 00283 IF( ANORM.GT.ZERO ) THEN 00284 RESULT( 1 ) = ( ( RESID / REAL( MAX(1,M,N) ) ) / ANORM ) / ULP 00285 ELSE 00286 RESULT( 1 ) = ZERO 00287 END IF 00288 * 00289 * Compute T*Z - Q'*B 00290 * 00291 CALL SGEMM( 'No Transpose', 'No transpose', N, P, P, ONE, T, LDB, 00292 $ Z, LDB, ZERO, BWK, LDB ) 00293 CALL SGEMM( 'Transpose', 'No transpose', N, P, N, -ONE, Q, LDA, 00294 $ B, LDB, ONE, BWK, LDB ) 00295 * 00296 * Compute norm( T*Z - Q'*B ) / ( MAX(P,N)*norm(A)*ULP ) . 00297 * 00298 RESID = SLANGE( '1', N, P, BWK, LDB, RWORK ) 00299 IF( BNORM.GT.ZERO ) THEN 00300 RESULT( 2 ) = ( ( RESID / REAL( MAX(1,P,N ) ) )/BNORM ) / ULP 00301 ELSE 00302 RESULT( 2 ) = ZERO 00303 END IF 00304 * 00305 * Compute I - Q'*Q 00306 * 00307 CALL SLASET( 'Full', N, N, ZERO, ONE, R, LDA ) 00308 CALL SSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDA, ONE, R, 00309 $ LDA ) 00310 * 00311 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00312 * 00313 RESID = SLANSY( '1', 'Upper', N, R, LDA, RWORK ) 00314 RESULT( 3 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP 00315 * 00316 * Compute I - Z'*Z 00317 * 00318 CALL SLASET( 'Full', P, P, ZERO, ONE, T, LDB ) 00319 CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, Z, LDB, ONE, T, 00320 $ LDB ) 00321 * 00322 * Compute norm( I - Z'*Z ) / ( P*ULP ) . 00323 * 00324 RESID = SLANSY( '1', 'Upper', P, T, LDB, RWORK ) 00325 RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP 00326 * 00327 RETURN 00328 * 00329 * End of SGQRTS 00330 * 00331 END