![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGRQTS 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 DGRQTS( M, P, N, 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, N, P 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00019 * $ BF( LDB, * ), BWK( LDB, * ), Q( LDA, * ), 00020 * $ R( LDA, * ), RESULT( 4 ), RWORK( * ), 00021 * $ T( LDB, * ), TAUA( * ), TAUB( * ), 00022 * $ WORK( LWORK ), Z( LDB, * ) 00023 * .. 00024 * 00025 * 00026 *> \par Purpose: 00027 * ============= 00028 *> 00029 *> \verbatim 00030 *> 00031 *> DGRQTS tests DGGRQF, which computes the GRQ factorization of an 00032 *> M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q. 00033 *> \endverbatim 00034 * 00035 * Arguments: 00036 * ========== 00037 * 00038 *> \param[in] M 00039 *> \verbatim 00040 *> M is INTEGER 00041 *> The number of rows of the matrix A. M >= 0. 00042 *> \endverbatim 00043 *> 00044 *> \param[in] P 00045 *> \verbatim 00046 *> P is INTEGER 00047 *> The number of rows of the matrix B. P >= 0. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] N 00051 *> \verbatim 00052 *> N is INTEGER 00053 *> The number of columns of the matrices A and B. N >= 0. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] A 00057 *> \verbatim 00058 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00059 *> The M-by-N matrix A. 00060 *> \endverbatim 00061 *> 00062 *> \param[out] AF 00063 *> \verbatim 00064 *> AF is DOUBLE PRECISION array, dimension (LDA,N) 00065 *> Details of the GRQ factorization of A and B, as returned 00066 *> by DGGRQF, see SGGRQF for further details. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] Q 00070 *> \verbatim 00071 *> Q is DOUBLE PRECISION array, dimension (LDA,N) 00072 *> The N-by-N orthogonal matrix Q. 00073 *> \endverbatim 00074 *> 00075 *> \param[out] R 00076 *> \verbatim 00077 *> R is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(M,N)) 00090 *> The scalar factors of the elementary reflectors, as returned 00091 *> by DGGQRC. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] B 00095 *> \verbatim 00096 *> B is DOUBLE PRECISION array, dimension (LDB,N) 00097 *> On entry, the P-by-N matrix A. 00098 *> \endverbatim 00099 *> 00100 *> \param[out] BF 00101 *> \verbatim 00102 *> BF is DOUBLE PRECISION array, dimension (LDB,N) 00103 *> Details of the GQR factorization of A and B, as returned 00104 *> by DGGRQF, see SGGRQF for further details. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] Z 00108 *> \verbatim 00109 *> Z is DOUBLE PRECISION array, dimension (LDB,P) 00110 *> The P-by-P orthogonal matrix Z. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] T 00114 *> \verbatim 00115 *> T is DOUBLE PRECISION array, dimension (LDB,max(P,N)) 00116 *> \endverbatim 00117 *> 00118 *> \param[out] BWK 00119 *> \verbatim 00120 *> BWK is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (min(P,N)) 00133 *> The scalar factors of the elementary reflectors, as returned 00134 *> by DGGRQF. 00135 *> \endverbatim 00136 *> 00137 *> \param[out] WORK 00138 *> \verbatim 00139 *> WORK is DOUBLE PRECISION 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(M,P,N)**2. 00146 *> \endverbatim 00147 *> 00148 *> \param[out] RWORK 00149 *> \verbatim 00150 *> RWORK is DOUBLE PRECISION array, dimension (M) 00151 *> \endverbatim 00152 *> 00153 *> \param[out] RESULT 00154 *> \verbatim 00155 *> RESULT is DOUBLE PRECISION array, dimension (4) 00156 *> The test ratios: 00157 *> RESULT(1) = norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP) 00158 *> RESULT(2) = norm( T*Q - Z'*B ) / (MAX(P,N)*norm(B)*ULP) 00159 *> RESULT(3) = norm( I - Q'*Q ) / ( N*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 double_eig 00174 * 00175 * ===================================================================== 00176 SUBROUTINE DGRQTS( M, P, N, 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, N, P 00186 * .. 00187 * .. Array Arguments .. 00188 DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), B( LDB, * ), 00189 $ BF( LDB, * ), BWK( LDB, * ), Q( LDA, * ), 00190 $ R( LDA, * ), RESULT( 4 ), RWORK( * ), 00191 $ T( LDB, * ), TAUA( * ), TAUB( * ), 00192 $ WORK( LWORK ), Z( LDB, * ) 00193 * .. 00194 * 00195 * ===================================================================== 00196 * 00197 * .. Parameters .. 00198 DOUBLE PRECISION ZERO, ONE 00199 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00200 DOUBLE PRECISION ROGUE 00201 PARAMETER ( ROGUE = -1.0D+10 ) 00202 * .. 00203 * .. Local Scalars .. 00204 INTEGER INFO 00205 DOUBLE PRECISION ANORM, BNORM, RESID, ULP, UNFL 00206 * .. 00207 * .. External Functions .. 00208 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 00209 EXTERNAL DLAMCH, DLANGE, DLANSY 00210 * .. 00211 * .. External Subroutines .. 00212 EXTERNAL DGEMM, DGGRQF, DLACPY, DLASET, DORGQR, DORGRQ, 00213 $ DSYRK 00214 * .. 00215 * .. Intrinsic Functions .. 00216 INTRINSIC DBLE, MAX, MIN 00217 * .. 00218 * .. Executable Statements .. 00219 * 00220 ULP = DLAMCH( 'Precision' ) 00221 UNFL = DLAMCH( 'Safe minimum' ) 00222 * 00223 * Copy the matrix A to the array AF. 00224 * 00225 CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA ) 00226 CALL DLACPY( 'Full', P, N, B, LDB, BF, LDB ) 00227 * 00228 ANORM = MAX( DLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) 00229 BNORM = MAX( DLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) 00230 * 00231 * Factorize the matrices A and B in the arrays AF and BF. 00232 * 00233 CALL DGGRQF( M, P, N, AF, LDA, TAUA, BF, LDB, TAUB, WORK, LWORK, 00234 $ INFO ) 00235 * 00236 * Generate the N-by-N matrix Q 00237 * 00238 CALL DLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA ) 00239 IF( M.LE.N ) THEN 00240 IF( M.GT.0 .AND. M.LT.N ) 00241 $ CALL DLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA ) 00242 IF( M.GT.1 ) 00243 $ CALL DLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA, 00244 $ Q( N-M+2, N-M+1 ), LDA ) 00245 ELSE 00246 IF( N.GT.1 ) 00247 $ CALL DLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA, 00248 $ Q( 2, 1 ), LDA ) 00249 END IF 00250 CALL DORGRQ( N, N, MIN( M, N ), Q, LDA, TAUA, WORK, LWORK, INFO ) 00251 * 00252 * Generate the P-by-P matrix Z 00253 * 00254 CALL DLASET( 'Full', P, P, ROGUE, ROGUE, Z, LDB ) 00255 IF( P.GT.1 ) 00256 $ CALL DLACPY( 'Lower', P-1, N, BF( 2, 1 ), LDB, Z( 2, 1 ), LDB ) 00257 CALL DORGQR( P, P, MIN( P, N ), Z, LDB, TAUB, WORK, LWORK, INFO ) 00258 * 00259 * Copy R 00260 * 00261 CALL DLASET( 'Full', M, N, ZERO, ZERO, R, LDA ) 00262 IF( M.LE.N ) THEN 00263 CALL DLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, R( 1, N-M+1 ), 00264 $ LDA ) 00265 ELSE 00266 CALL DLACPY( 'Full', M-N, N, AF, LDA, R, LDA ) 00267 CALL DLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, R( M-N+1, 1 ), 00268 $ LDA ) 00269 END IF 00270 * 00271 * Copy T 00272 * 00273 CALL DLASET( 'Full', P, N, ZERO, ZERO, T, LDB ) 00274 CALL DLACPY( 'Upper', P, N, BF, LDB, T, LDB ) 00275 * 00276 * Compute R - A*Q' 00277 * 00278 CALL DGEMM( 'No transpose', 'Transpose', M, N, N, -ONE, A, LDA, Q, 00279 $ LDA, ONE, R, LDA ) 00280 * 00281 * Compute norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP ) . 00282 * 00283 RESID = DLANGE( '1', M, N, R, LDA, RWORK ) 00284 IF( ANORM.GT.ZERO ) THEN 00285 RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / 00286 $ ULP 00287 ELSE 00288 RESULT( 1 ) = ZERO 00289 END IF 00290 * 00291 * Compute T*Q - Z'*B 00292 * 00293 CALL DGEMM( 'Transpose', 'No transpose', P, N, P, ONE, Z, LDB, B, 00294 $ LDB, ZERO, BWK, LDB ) 00295 CALL DGEMM( 'No transpose', 'No transpose', P, N, N, ONE, T, LDB, 00296 $ Q, LDA, -ONE, BWK, LDB ) 00297 * 00298 * Compute norm( T*Q - Z'*B ) / ( MAX(P,N)*norm(A)*ULP ) . 00299 * 00300 RESID = DLANGE( '1', P, N, BWK, LDB, RWORK ) 00301 IF( BNORM.GT.ZERO ) THEN 00302 RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, M ) ) ) / BNORM ) / 00303 $ ULP 00304 ELSE 00305 RESULT( 2 ) = ZERO 00306 END IF 00307 * 00308 * Compute I - Q*Q' 00309 * 00310 CALL DLASET( 'Full', N, N, ZERO, ONE, R, LDA ) 00311 CALL DSYRK( 'Upper', 'No Transpose', N, N, -ONE, Q, LDA, ONE, R, 00312 $ LDA ) 00313 * 00314 * Compute norm( I - Q'*Q ) / ( N * ULP ) . 00315 * 00316 RESID = DLANSY( '1', 'Upper', N, R, LDA, RWORK ) 00317 RESULT( 3 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP 00318 * 00319 * Compute I - Z'*Z 00320 * 00321 CALL DLASET( 'Full', P, P, ZERO, ONE, T, LDB ) 00322 CALL DSYRK( 'Upper', 'Transpose', P, P, -ONE, Z, LDB, ONE, T, 00323 $ LDB ) 00324 * 00325 * Compute norm( I - Z'*Z ) / ( P*ULP ) . 00326 * 00327 RESID = DLANSY( '1', 'Upper', P, T, LDB, RWORK ) 00328 RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP 00329 * 00330 RETURN 00331 * 00332 * End of DGRQTS 00333 * 00334 END