![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SQRT03 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 SQRT03( M, N, K, AF, C, CC, Q, LDA, TAU, WORK, LWORK, 00012 * RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * INTEGER K, LDA, LWORK, M, N 00016 * .. 00017 * .. Array Arguments .. 00018 * REAL AF( LDA, * ), C( LDA, * ), CC( LDA, * ), 00019 * $ Q( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), 00020 * $ WORK( LWORK ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SQRT03 tests SORMQR, which computes Q*C, Q'*C, C*Q or C*Q'. 00030 *> 00031 *> SQRT03 compares the results of a call to SORMQR with the results of 00032 *> forming Q explicitly by a call to SORGQR and then performing matrix 00033 *> multiplication by a call to SGEMM. 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] M 00040 *> \verbatim 00041 *> M is INTEGER 00042 *> The order of the orthogonal matrix Q. M >= 0. 00043 *> \endverbatim 00044 *> 00045 *> \param[in] N 00046 *> \verbatim 00047 *> N is INTEGER 00048 *> The number of rows or columns of the matrix C; C is m-by-n if 00049 *> Q is applied from the left, or n-by-m if Q is applied from 00050 *> the right. N >= 0. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] K 00054 *> \verbatim 00055 *> K is INTEGER 00056 *> The number of elementary reflectors whose product defines the 00057 *> orthogonal matrix Q. M >= K >= 0. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] AF 00061 *> \verbatim 00062 *> AF is REAL array, dimension (LDA,N) 00063 *> Details of the QR factorization of an m-by-n matrix, as 00064 *> returnedby SGEQRF. See SGEQRF for further details. 00065 *> \endverbatim 00066 *> 00067 *> \param[out] C 00068 *> \verbatim 00069 *> C is REAL array, dimension (LDA,N) 00070 *> \endverbatim 00071 *> 00072 *> \param[out] CC 00073 *> \verbatim 00074 *> CC is REAL array, dimension (LDA,N) 00075 *> \endverbatim 00076 *> 00077 *> \param[out] Q 00078 *> \verbatim 00079 *> Q is REAL array, dimension (LDA,M) 00080 *> \endverbatim 00081 *> 00082 *> \param[in] LDA 00083 *> \verbatim 00084 *> LDA is INTEGER 00085 *> The leading dimension of the arrays AF, C, CC, and Q. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] TAU 00089 *> \verbatim 00090 *> TAU is REAL array, dimension (min(M,N)) 00091 *> The scalar factors of the elementary reflectors corresponding 00092 *> to the QR factorization in AF. 00093 *> \endverbatim 00094 *> 00095 *> \param[out] WORK 00096 *> \verbatim 00097 *> WORK is REAL array, dimension (LWORK) 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LWORK 00101 *> \verbatim 00102 *> LWORK is INTEGER 00103 *> The length of WORK. LWORK must be at least M, and should be 00104 *> M*NB, where NB is the blocksize for this environment. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] RWORK 00108 *> \verbatim 00109 *> RWORK is REAL array, dimension (M) 00110 *> \endverbatim 00111 *> 00112 *> \param[out] RESULT 00113 *> \verbatim 00114 *> RESULT is REAL array, dimension (4) 00115 *> The test ratios compare two techniques for multiplying a 00116 *> random matrix C by an m-by-m orthogonal matrix Q. 00117 *> RESULT(1) = norm( Q*C - Q*C ) / ( M * norm(C) * EPS ) 00118 *> RESULT(2) = norm( C*Q - C*Q ) / ( M * norm(C) * EPS ) 00119 *> RESULT(3) = norm( Q'*C - Q'*C )/ ( M * norm(C) * EPS ) 00120 *> RESULT(4) = norm( C*Q' - C*Q' )/ ( M * norm(C) * EPS ) 00121 *> \endverbatim 00122 * 00123 * Authors: 00124 * ======== 00125 * 00126 *> \author Univ. of Tennessee 00127 *> \author Univ. of California Berkeley 00128 *> \author Univ. of Colorado Denver 00129 *> \author NAG Ltd. 00130 * 00131 *> \date November 2011 00132 * 00133 *> \ingroup single_lin 00134 * 00135 * ===================================================================== 00136 SUBROUTINE SQRT03( M, N, K, AF, C, CC, Q, LDA, TAU, WORK, LWORK, 00137 $ RWORK, RESULT ) 00138 * 00139 * -- LAPACK test routine (version 3.4.0) -- 00140 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00141 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00142 * November 2011 00143 * 00144 * .. Scalar Arguments .. 00145 INTEGER K, LDA, LWORK, M, N 00146 * .. 00147 * .. Array Arguments .. 00148 REAL AF( LDA, * ), C( LDA, * ), CC( LDA, * ), 00149 $ Q( LDA, * ), RESULT( * ), RWORK( * ), TAU( * ), 00150 $ WORK( LWORK ) 00151 * .. 00152 * 00153 * ===================================================================== 00154 * 00155 * .. Parameters .. 00156 REAL ONE 00157 PARAMETER ( ONE = 1.0E0 ) 00158 REAL ROGUE 00159 PARAMETER ( ROGUE = -1.0E+10 ) 00160 * .. 00161 * .. Local Scalars .. 00162 CHARACTER SIDE, TRANS 00163 INTEGER INFO, ISIDE, ITRANS, J, MC, NC 00164 REAL CNORM, EPS, RESID 00165 * .. 00166 * .. External Functions .. 00167 LOGICAL LSAME 00168 REAL SLAMCH, SLANGE 00169 EXTERNAL LSAME, SLAMCH, SLANGE 00170 * .. 00171 * .. External Subroutines .. 00172 EXTERNAL SGEMM, SLACPY, SLARNV, SLASET, SORGQR, SORMQR 00173 * .. 00174 * .. Local Arrays .. 00175 INTEGER ISEED( 4 ) 00176 * .. 00177 * .. Intrinsic Functions .. 00178 INTRINSIC MAX, REAL 00179 * .. 00180 * .. Scalars in Common .. 00181 CHARACTER*32 SRNAMT 00182 * .. 00183 * .. Common blocks .. 00184 COMMON / SRNAMC / SRNAMT 00185 * .. 00186 * .. Data statements .. 00187 DATA ISEED / 1988, 1989, 1990, 1991 / 00188 * .. 00189 * .. Executable Statements .. 00190 * 00191 EPS = SLAMCH( 'Epsilon' ) 00192 * 00193 * Copy the first k columns of the factorization to the array Q 00194 * 00195 CALL SLASET( 'Full', M, M, ROGUE, ROGUE, Q, LDA ) 00196 CALL SLACPY( 'Lower', M-1, K, AF( 2, 1 ), LDA, Q( 2, 1 ), LDA ) 00197 * 00198 * Generate the m-by-m matrix Q 00199 * 00200 SRNAMT = 'SORGQR' 00201 CALL SORGQR( M, M, K, Q, LDA, TAU, WORK, LWORK, INFO ) 00202 * 00203 DO 30 ISIDE = 1, 2 00204 IF( ISIDE.EQ.1 ) THEN 00205 SIDE = 'L' 00206 MC = M 00207 NC = N 00208 ELSE 00209 SIDE = 'R' 00210 MC = N 00211 NC = M 00212 END IF 00213 * 00214 * Generate MC by NC matrix C 00215 * 00216 DO 10 J = 1, NC 00217 CALL SLARNV( 2, ISEED, MC, C( 1, J ) ) 00218 10 CONTINUE 00219 CNORM = SLANGE( '1', MC, NC, C, LDA, RWORK ) 00220 IF( CNORM.EQ.0.0 ) 00221 $ CNORM = ONE 00222 * 00223 DO 20 ITRANS = 1, 2 00224 IF( ITRANS.EQ.1 ) THEN 00225 TRANS = 'N' 00226 ELSE 00227 TRANS = 'T' 00228 END IF 00229 * 00230 * Copy C 00231 * 00232 CALL SLACPY( 'Full', MC, NC, C, LDA, CC, LDA ) 00233 * 00234 * Apply Q or Q' to C 00235 * 00236 SRNAMT = 'SORMQR' 00237 CALL SORMQR( SIDE, TRANS, MC, NC, K, AF, LDA, TAU, CC, LDA, 00238 $ WORK, LWORK, INFO ) 00239 * 00240 * Form explicit product and subtract 00241 * 00242 IF( LSAME( SIDE, 'L' ) ) THEN 00243 CALL SGEMM( TRANS, 'No transpose', MC, NC, MC, -ONE, Q, 00244 $ LDA, C, LDA, ONE, CC, LDA ) 00245 ELSE 00246 CALL SGEMM( 'No transpose', TRANS, MC, NC, NC, -ONE, C, 00247 $ LDA, Q, LDA, ONE, CC, LDA ) 00248 END IF 00249 * 00250 * Compute error in the difference 00251 * 00252 RESID = SLANGE( '1', MC, NC, CC, LDA, RWORK ) 00253 RESULT( ( ISIDE-1 )*2+ITRANS ) = RESID / 00254 $ ( REAL( MAX( 1, M ) )*CNORM*EPS ) 00255 * 00256 20 CONTINUE 00257 30 CONTINUE 00258 * 00259 RETURN 00260 * 00261 * End of SQRT03 00262 * 00263 END